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Abstract 


Simulations  of  muon  interactions  with  high  Z  material  using  two  different  muon 
energies,  100  MeV  and  1  GeV,  were  performed  on  Eve  different  materials  of  various 
atomic  numbers  yielding  average  neutron  production  rates  that  range  from  2.3  ± 
0.01  in  enriched  uranium  to  negligible  amounts  in  aluminum  when  exposed  to  the 
100  MeV  energy  muons.  As  the  muon  energy  was  increased  to  1  GeV,  neutron  yields 
shrank  to  negligible  levels.  Little  difference  was  found  in  neutron  yield  produced  in 
non-hssile  material. 

Experimental  data  was  collected  by  exposing  a  15  cm  thick  block  of  iron,  and  5 
and  15  cm  thick  blocks  of  lead  to  the  natural  atmospheric  muon  flux.  The  incident 
muon  energy  distribution  was  found  to  have  a  mode  of  180  MeV  and  a  mean  of 
520  MeV.  Probability  distributions  were  constructed  for  the  neutron  yields  of  each 
incident  muon  and  no  difference  was  found  in  the  various  distributions.  The  average 
muon  induced  neutron  yield  was  also  calculated  and  found  to  be  3.4  ±  0.1  for  a  15 
cm  thick  block  of  iron,  2.8  ±  0.1  for  a  5  cm  thick  block  of  lead,  and  2.2  ±  0.1  for  a 
15  cm  thick  block  of  lead. 
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ANALYSIS  OF  MUON  INDUCED  NEUTRONS  IN  DETECTING  HIGH  Z 


NUCLEAR  MATERIALS 

I.  Introduction 


1.1  Objective 

Currently  a  need  exists  to  detect  shielded  special  nuclear  material  using  passive 
interrogation  techniques.  Muon  imaging  holds  the  potential  to  satisfy  this  require¬ 
ment,  but  current  implementations  are  hindered  by  long  detection  times.  A  proposed 
solution  is  to  utilize  neutrons  produced  by  muon  interactions  with  the  fissile  material 
to  supplement  current  muon  imaging  techniques.  It  is  predicted  that  muon  induced 
fission  of  fissile  nuclear  materials  will  produce  enough  detectiblc  neutrons  through 
interactions  within  the  target  that  resulting  signals  can  be  utilized  as  an  additional 
detection  method  when  neutrons  are  detected  in  coincidence  with  incoming  muons. 

The  primary  objective  of  this  research  is  to  examine  the  neutron  production  rate 
of  various  materials  when  exposed  to  atmospheric  muons.  As  not  all  material  con¬ 
tained  within  a  nuclear  weapon  is  fissile,  it  is  important  to  understand  the  neutron 
production  capabilities  of  a  wide  range  of  atomic  number  elements.  It  is  anticipated 
that  the  actinide  elements,  especially  fissile  actinides  which  have  a  low  activation 
energy  for  fission,  will  produce  a  higher  number  of  neutrons  through  muon  induced 
fission  events,  while  neutron  production  will  decrease  for  the  lower  Z  materials.  This 
research  will  focus  on  the  muon  induced  neutron  production  from  both  iron  and  lead. 
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1.2  Motivation 


Accurate  detection  of  Special  Nuclear  Material  (SNM)  is  of  paramount  importance 
to  issues  of  national  security  and  international  treaty  verification.  With  the  rise  of 
non-state  actor  terrorist  organizations,  a  fear  of  a  nuclear  weapon  falling  into  the 
hands  of  one  of  these  organizations  and  being  smuggled  into  the  United  States  has 
become  more  viable.  This  fear  has  prompted  the  Department  of  Homeland  Security 
(DHS)  and  Department  of  Defense  (DoD)  to  work  with  the  Department  of  Energy 
(DoE)  laboratories  to  create  detection  portals  for  use  at  locations  such  as  border 
crossings  and  harbors.  These  portals  utilize  detectors  such  as  thallium  doped  sodium 
iodide  (Nal(Tl))  scintillation  detectors  and  High  Purity  Germanium  (HPGe)  detectors 
to  indirectly  inspect  incoming  vehicles  and/or  shipping  containers  for  any  type  of 
potential  nuclear  material  or  Radiological  Dispersion  Device  (RDD).  However,  these 
detection  methods  can  be  defeated  with  ample  shielding  of  the  smuggled  material. 
If  surrounded  by  enough  material  at  a  high  atomic  number,  the  gamma  rays  can  be 
sufficiently  attenuated  to  undetectable  levels. 

Moreover,  as  the  nuclear  arsenals  of  the  major  nuclear  weapon  states  continue  to 
decrease,  the  need  arises  for  more  accurate  verification  methods  of  the  material  con¬ 
tained  within  these  weapons.  Visual  inspections  and  active  interrogation  techniques 
are  prohibited,  due  to  the  possibility  of  revealing  sensitive  weapon  design  information. 
Additionally,  thick  outer  casing  material  can  have  the  same  shielding  problem  in  the 
portal  scenario.  Since  active  examination  of  internal  components  are  not  allowed,  an 
alternate  method  of  imaging  and  detecting  special  nuclear  material  is  required.  Muon 
imaging  has  been  suggested  as  one  solution  to  this  problem  [1], 

Muon  tomography  utilizes  negatively  charged  muons  created  by  cosmic  particle 
interactions  in  the  upper  atmosphere.  These  particles  are  207  times  more  massive 
than  an  electron  and  have  a  mean  lifetime  of  2.2  /xs  [2],  Because  of  their  capability  to 
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penetrate  through  dense  materials,  muons  have  been  used  to  image  major  geograph¬ 
ical  landmarks.  In  one  instance  they  were  utilized  to  examine  the  magma  buildup 
in  an  active  volcano  located  in  Japan  [3].  In  another,  researchers  capitalized  on  the 
muon’s  penetration  capability  to  examine  the  pyramids  for  hidden  chambers  [4].  A 
similar  technique  could  be  applied  to  weapons  inspection  procedures  to  probe  the 
unobservable  physics  package  for  dense  nuclear  fuel  even  in  the  presence  of  dense 
shielding  material. 

When  muon  tomography  alone  is  used,  the  time  period  required  for  the  necessary 
resolution  of  meaningful  information  is  on  the  order  of  tens  of  hours  [5] .  However,  it 
may  be  possible  to  decrease  the  detection  time  of  nuclear  material  if  muon  tomography 
is  combined  with  neutron  detection  from  muon  induced  fission.  In  this  process,  a 
muon  may  interact  with  a  heavy  nucleus  and  replace  an  electron  in  the  atom.  As 
the  muon  deexcites  and  is  captured  by  the  heavy  nucleus,  it  has  a  chance  to  transfer 
its  energy  to  the  nucleus  which  can  be  enough  to  overcome  the  fission  activation 
barrier  if  the  fission  barrier  is  low.  These  fission  events  produce  excess  neutrons  in 
the  system  which  can  be  detected  and  utilized  in  further  analysis  of  the  system  under 
examination. 
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II.  Theory 


As  stated  previously,  muons  are  subatomic  particles  known  as  leptons  with  a 
mass  207  times  that  of  the  electron  and  can  have  either  a  positive  or  negative  charge. 
They  are  primarily  produced  in  the  upper  atmosphere  through  collisions  of  cosmic 
rays,  predominately  very  energetic  protons,  and  the  atmospheric  molecules.  As  these 
high  energy  protons  collide  with  other  particles,  pions  are  produced  which  quickly 
decay  into  both  positive  and  negative  muons  [6].  These  are  much  longer  lived,  with 
a  mean  lifetime  of  2.2  /rs.  Since  they  are  created  at  such  high  energies  and,  thus 
having  relativistic  velocities,  these  particles  can  easily  reach  the  surface  of  the  earth 
to  interact  with  matter  there.  An  example  of  this  creation  process  can  be  seen  in 
Figure  1.  These  muons  will  reach  the  surface  at  a  rate  of  approximately  10,000 
muons  /  min  /  m2 . 


e+ 


Figure  1.  High  energy  collisions  of  cosmic  ray  protons  in  the  upper  atmosphere  will 
produce  pions  that  quickly  decay  into  muons  among  other  perticles. 
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2.1  Muon  Imaging 


Because  of  the  nearly  constant  muon  flux,  these  naturally  occurring  muons  have 
been  used  to  image  various  objects  since  the  early  1970s.  At  that  time  Alvarez  et 
al.  capitalized  on  the  muons  ability  to  penetrate  through  large  amounts  of  material 
by  imaging  Cepherens  pyramid  [4],  In  this  particular  experiment,  the  inspecting 
team  placed  detectors  beneath  the  desired  location  to  be  inspected  and  gathered 
data  for  several  months.  Once  the  exposure  was  completed,  the  detected  muon  flux 
was  compared  to  the  known  background  levels.  By  accounting  for  the  anticipated 
attenuation  of  the  muons  in  the  density  of  material  overhead,  the  internal  structure 
of  the  pyramid  could  be  mapped  out  to  rule  out  the  possibility  of  a  hidden  chamber 
within  Cepherens  pyramid  which  had  been  previously  found  in  others. 

In  another  similar  experiment,  Tanaka,  et  al.  placed  one  4000  cm2  AgBr  emulsion 
cloud  chamber  located  underground  approximately  1  km  away  from  the  Asama  vol¬ 
cano  in  Japan  to  gather  density  data  [3].  To  confirm  the  accuracy  of  their  findings, 
the  team  then  imaged  the  Usu  lava  dome  and  compared  to  other  known  density  mea¬ 
surements.  These  findings  can  be  seen  in  Figure  2.  This  same  technique  has  also  been 
suggested  for  application  in  imaging  damaged  nuclear  reactors  such  as  Fukushima, 
that  it  would  be  unsafe  to  approach  and  inspect  at  a  close  distance  [7]  [8]. 

In  both  of  these  illustrations,  the  muon  imaging  was  conducted  over  very  long 
periods  of  time  to  observe  a  sufficient  decrease  in  the  muon  flux  due  to  attenuation  in 
the  dense  material.  Two  other  forms  of  muon  imaging  also  exist.  Muon  tomography 
works  by  measuring  the  amount  of  scattered  muons.  Two  detectors  are  placed  on 
either  side  of  the  object,  one  to  measure  the  incoming  angle  of  the  muon,  and  another 
to  measure  the  exiting  direction.  Since  muons  are  much  more  likely  to  scatter  at  a 
greater  angle  when  passing  through  dense  material,  an  idea  of  the  location  of  the 
material  in  between  the  two  detectors  can  be  gathered.  In  2007  an  experiment  was 
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Muon  Detector  0  150  aoom 


Figure  2.  Cosmic  muon  imaging  results  provided  the  density  measurements  for  the  lava 
dome.  Alternate  density  measurements  were  available  and  confirmed  the  accuracy  of 
the  muon  imaging  technique.  [3]  This  image  has  been  reproduced  with  the  permission 
of  the  American  Journal  of  Science. 

conducted  at  CERN  to  examine  this  very  technique.  In  the  experiment,  two  lead 
blocks  were  placed  between  gas  filled  drift  tubes  in  the  setup  seen  in  Figure  3  [9]. 
The  results  of  the  CERN  experiment  can  be  seen  in  Figure  4. 

Further  studies  have  been  conducted  using  this  same  technique  and  have  produced 
promising  results.  The  DOE  has  even  worked  with  Decision  Sciences  International 
Corporation  to  develop  portal  systems  which  have  been  shown  to  clear  a  cargo  con¬ 
tainer  in  approximately  30  seconds  [10]. 

Another  method  of  muon  imaging  is  known  as  telescopic  mode.  In  this  method, 
it  is  only  necessary  to  detect  incoming  muons.  Of  all  the  muons  that  enter  the 
material,  some  will  be  moderated  to  a  point  that  they  will  stop  and  be  captured 
within  the  target  material.  At  this  point,  secondary  signals  can  be  given  off  which 
can  be  detected  to  confirm  the  presence  of  various  types  of  material.  To  confirm  the 
presence  of  nuclear  material  specifically,  a  secondary  signal  must  be  acquired  through 
either  gamma  or  neutron  production  via  muon  induced  fission  events.  Several  studies 
have  been  conducted  to  examine  the  feasibility  of  this  particular  method  and  have 
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Figure  3.  Experimental  setup  utilized  in  the  muon  tomography  experiment  run  by 
CERN  when  inspecting  the  location  of  two  lead  blocks.  [9]  Copyright  (2007)  IEEE 


Figure  4.  Experimental  results  demonstrate  muon  tomography’s  capability  to  detect 
and  determine  location  of  lead  blocks.  [9]  Copyright  (2007)  IEEE 
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shown  much  promise  [6].  It  is  this  type  of  muon  imaging,  telescopic  mode,  that  will 
be  the  focus  of  this  experiment. 

It  should  be  noted  that  muon  tomography  and  telescopic  mode  are  not  mutu¬ 
ally  exclusive.  In  fact,  they  may  be  complimentary  to  one  another  based  upon  the 
muon  energy  dependence  of  the  phenomena.  High  energy  muons  primarily  interact 
at  a  localized  point  via  scattering  events  which  can  be  utilized  through  scattering 
tomography.  Lower  energy  muons  have  a  much  higher  linear  energy  transfer  making 
moderation  within  the  material  much  more  probable  and  increasing  the  effectiveness 
of  neutron  production  and  gamma  emission.  By  using  both  methods  simultaneously, 
more  information  can  be  gathered  about  the  target  material  than  by  using  one  method 
exclusively. 

2.2  Muon  Induced  Fission 

Fission  occurs  when  a  nucleus  is  imparted  with  enough  energy  to  overcome  the 
strong  nuclear  force  holding  its  nucleons  together.  When  this  occurs,  the  nucleus  will 
split  into  two  daughter  nuclei  and  several  neutrons,  and  the  masses  of  the  individual 
components  will  sum  to  less  than  that  of  the  original.  This  change  in  mass  is  released 
as  energy  according  to  the  equation,  E  =  me2.  The  fission  phenomenon  is  well 
understood  by  the  scientific  community  and  has  been  utilized  in  both  the  creation  of 
weapons  and  power  plants  for  energy  production. 

Muon  Induced  Fission  (//IF)  was  originally  proposed  by  Wheeler  in  1948  [11]. 
Since  that  time  it  has  been  explored  in  much  more  depth  and  has  proven  useful  in 
studying  nuclear  energy  dissipation  and  fission  dynamics  [12]  [13].  //IF  occurs  by 
two  processes,  by  electromagnetic  radiationless  transition  causing  prompt  fission  or 
by  weak  capture  reactions  causing  delayed  fission.  In  both  cases,  it  is  important  to 
note  the  time  scales  of  muon  capture  occur  on  timescales  that  are  small  compared  to 


the  mean  lifetime  of  the  muon  itself  [14]  [15].  This  allows  the  muons  to  be  captured 
by  the  material  and  transfer  the  energy  to  a  nucleus.  If  the  energy  transfer  is  large 
enough,  the  fission  barrier  can  be  overcome  and  the  nucleus  will  split.  Alternatively,  if 
the  energy  transfer  is  not  sufficient  for  fission,  the  muon  can  be  captured  by  a  proton 
within  the  nucleus,  converting  the  atom  to  another  element. 

To  initiate  the  fission  process,  a  negatively  charged  muon  must  first  interact  with 
the  material.  As  the  muon  enters  the  fissionable  material  and  begins  to  interact 
through  collisions,  it  rapidly  loses  a  majority  of  its  energy  by  ionizing  the  surrounding 
material  within  10-9  to  lCT10s  [15].  This  process  may  continue  until  enough  time  has 
passed  and  the  muon  eventually  reaches  the  end  of  its  lifetime.  When  this  occurs,  the 
negative  muon  will  decay  into  a  muon  neutrino,  an  electron,  and  an  electron  neutrino, 
and  the  reaction  can  be  seen  in  Equation  1. 

h  — >■  ^  +  e  +  Pe  (1) 

Because  the  muon  decay  process  is  slow  compared  to  the  processes  under  con¬ 
sideration,  its  effects  are  negligible.  An  alternative  to  this  scenario  occurs  with  the 
capture  of  a  muon.  In  such  an  event,  the  muon  will  replace  an  electron  in  one  of  the 
outer  electron  orbitals  forming  an  excited  muonic  atom  [15].  From  this  position,  /df 
in  both  prompt  and  delayed  processes  may  occur. 

Prompt  Muon  Induced  Fission. 

In  the  case  of  ^IF,  prompt  fission  occurs  through  the  energy  transfer  of  the  muons’ 
transitions  within  the  atomic  energy  levels.  Because  a  muon  acts  as  a  heavy  electron, 
its  behavior  can  be  thought  of  in  an  analogous  manner.  Once  captured  the  muon  will 
reside  in  the  outer  orbitals  of  the  atom  until  it  deexcites  to  one  of  the  lower  energy 
levels.  The  2p—ls  and  3d  —  Is  transitions  are  on  the  order  of  the  fission  barrier  in  the 
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actinide  elements  [16]  [15].  In  these  situations,  the  energy  can  be  transferred  directly 
to  the  nucleus  rather  than  expelled  from  the  atom  via  electromagnetic  radiation  in 
a  process  known  as  inverse  internal  conversion  [15].  After  such  an  event  has  taken 
place,  the  nucleus  will  be  imparted  with  enough  energy  to  exceed  the  fission  barrier. 
The  nucleus  will  split  and,  most  often,  the  muon  will  remain  attached  to  the  larger 
fission  fragment.  However,  it  has  been  observed  that  under  certain  circumstances  the 
muon  will  attach  to  the  lighter  fission  fragment.  In  such  case  the  muon  attachment 
rates  to  the  lighter  daughter  product  can  be  examined  as  an  analysis  tool  to  study 
the  prompt  fission  dynamics  [15]. 

Delayed  Muon  Induced  Fission. 

Delayed  /xIF  requires  the  muon  to  have  been  captured  within  the  atomic  orbitals 
and  to  have  decayed  into  the  ground  state.  In  this  configuration  the  muon  will  spend 
a  significant  portion  of  its  time  residing  inside  the  nucleus  due  to  its  excessive  mass. 
At  that  location,  the  muon  can  be  captured  by  a  proton  resulting  in  a  neutron  and 
a  muon  neutrino.  This  process  yields  an  average  nucleus  excitation  energy  of  15-20 
MeV  which  is  much  greater  than  the  actinide  fission  barrier  of  5-6  MeV  [15].  Such 
a  large  energy  transfer  has  the  ability  to  allow  for  secondary  or  even  higher  fission 
events  [16] .  It  should  be  noted  that  the  mean  lifetime  of  these  events  is  based  upon  the 
weak  decay  process  on  the  order  of  10~8s  [15].  As  these  timescales  are  two  orders  of 
magnitude  less  than  the  muon  lifetime,  muon  decay  can  be  neglected  in  describing  the 
fission  processes.  Measurements  have  been  made  of  the  daughter  mass  distributions, 
delayed  fission  probabilities,  and  prompt  to  delayed  fission  ratios  [16]. 
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2.3  Muon  Catalyzed  Fusion 


Fusion  occurs  when  two  light  nuclei  combine  to  form  one  larger  nucleus  with  a 
mass  less  than  the  sum  of  the  original  components.  The  difference  in  mass  is  released 
as  energy  and  can  be  calculated  using  Einstein’s  equation,  E  =  me2.  In  practice, 
fusion  is  difficult  to  achieve  due  to  the  highly  repulsive  Coulombic  force  of  both 
positively  charged  nuclei.  Typically,  fusion  reactions  require  high  temperatures,  laser 
stimulation,  and  magnetic  confinement. 

However,  fusion  of  heavy  hydrogen  isotopes  can  be  achieved  without  the  need 
for  these  high  temperatures  or  laser  stimulation  by  way  of  Muon  Catalyzed  Fusion 
(/iCF).  In  this  process,  a  heavy  negatively  charged  muon  replaces  the  electron.  The 
muon  resides  between  two  hydrogen  isotopes  and  masks  the  coulomb  barrier  thereby 
reducing  the  atom’s  radius  by  over  200  times  that  of  a  normal  hydrogen  molecule  from 
approximately  10-8  cm  to  5.1xl0-11cm.  [17]  [18]  The  reduction  in  size  enables  higher 
rates  of  quantum  tunneling  of  one  nucleus  through  the  coloumb  barrier  of  the  other 
resulting  in  a  fusion  event.  This  technique  was  first  proposed  by  Andrei  Sakharov 
and  F.C.  Frank  in  1947  and  has  since  been  proven  experimentally  [19].  The  most 
common  fusion  reactions  are  those  of  Deuterium-Deuterium  (DD)  and  Deuterium- 
Tritium  (DT),  with  the  cross  section  for  the  DT  molecule  approximately  100  times 
larger  than  the  DD  molecule  [20]  [21]  .  The  entire  cycle  for  deuterium-deuterium 
fusion  can  be  seen  in  Figure  5. 

2.4  Muon  Spallation 

A  third  means  of  muon  neutron  production  is  via  muon  spallation.  In  this  method, 
the  incoming  muon  exchanges  a  virtual  photon  with  the  interaction  nucleus.  This  in¬ 
teraction  can  provide  the  nucleus  with  enough  energy  to  emit  one  or  several  neutrons. 
The  Feynman  diagram  in  Figure  6(a)  demonstrates  this  process  [22],  It  is  the  main 
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Figure  5.  The  cycle  of  muon  catalyzed  fusion  for  a  deuterium-deuterium  reaction. 
Reproduced  from  Dr.  Van  Dyk’s  dissertation.  [1] 


contribution  to  background  noise  in  deep  underground  experiments  searching  for  dark 
matter  and  neutrinoless  double  beta  decays  [23]. 


2.5  Neutron  Detectors 

In  order  to  properly  associate  neutron  events  with  an  incident  muon,  it  is  necessary 
to  construct  a  system  in  which  the  produced  neutrons  can  be  counted  in  coincidence 
with  a  muon  event.  To  accomplish  this  objective,  a  nearly  47t  solid  angle  neutron 
detection  system  should  be  implemented  to  produce  maximum  efficiency  of  neutron 
detection.  It  should  have  the  capability  to  discriminate  between  neutrons  and  gamma 
events  and  also  have  a  high  probability  of  interaction.  All  three  of  these  objectives 
can  be  accomplished  using  neutron/gamma  discriminating  liquid  scintillators.  De¬ 
tectors  like  the  BC-501A/EJ-301  and  BC-519/EJ-309  can  be  formed  into  the  desired 
configuration  in  order  to  gain  as  close  to  47t  coverage  as  possible  while  simultaneously 
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Figure  6.  (a) Feynman  diagram  depicting  neutron  production  through  the  exchange 

of  a  virtual  photon  between  a  muon  and  a  nucleus.  (b)Feynman  diagram  displaying 
neutron  production  through  photon  absorption.  Reproduced  from  Lawrence  Livermore 
National  Laboratory’s  report  on  Neutron  Production  by  Muon  Spallation.  [22] 


providing  the  capability  to  discriminate  between  neutron  events  and  gamma  events 
by  way  of  the  detection  method.  Neutron  events  are  detected  through  recoil  proton 
collisions  in  the  hydrogen  present  within  the  scintillation  material.  These  recoil  pro¬ 
tons  have  a  decay  time  which  is  much  longer  than  a  gamma  interaction  within  the 
same  material.  These  decay  times  are  associated  with  a  Pulse  Shape  Discrimination 
(PSD)  value  by  integrating  the  charge  collected  in  two  different  time  characteristics 
of  the  pulse,  Q short  and  Qiong ■  Both  of  these  charge  integrations  begin  at  the  trigger 
point,  and  can  then  be  used  to  calculate  the  PSD  value  by  using  Equation  2. 


PSD  = 


Qlcmg  Q  short 


Ql 


(2) 


ong 


By  plotting  a  2D  Histogram  of  the  PSD  values  calculated  by  electronic  digitizers 
after  the  exposure  of  the  detectors  to  a  neutron  source,  such  as  a  plutonium  beryllium 
source,  the  difference  in  PSD  values  for  neutrons  and  gammas  can  be  utilized  to 
distinguish  between  neutron  and  gamma  detections.  An  example  plot  can  be  seen 
in  Figure  7.  Once  a  similar  plot  has  been  produced,  cut  lines  can  be  determined 
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to  filter  out  the  extraneous  information  produced  from  gamma  detection  to  ensure 
that  only  neutron  events  are  counted  in  coincidence  with  incident  muons  within  the 
detection  system.  PSD  firmware,  such  as  that  produced  by  Caen  S.p.A,  is  utilized 
to  accomplish  all  these  tasks.  It  should  also  be  noted  that  there  is  an  overlap  of  the 
PSD  values  for  neutrons  and  gammas  at  low  energies.  This  overlap  is  not  entirely 
removed  when  only  a  single  cut  criterion  is  applied. 
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Figure  7.  A  2D  Histogram  of  PSD  values  shows  the  difference  between  neutron  and 
gamma  energy  depositions  within  the  liquid  scintillators.  [1] 

Detection  efficiency  also  must  be  considered  when  arranging  the  liquid  scintilla¬ 
tors.  The  arrangement  must  maximize  the  absolute  detector  efficiency  which  considers 
the  geometry  of  the  setup  as  well  as  the  intrinsic  efficiency  of  the  detector  itself.  The 
geometry  component  is  determined  by  the  solid  angle  of  the  source  exposed  to  the 
detector  when  compared  to  a  solid  angle  of  47t.  If  a  right  cylindrical  detector  is  con¬ 
sidered  and  a  point  source  is  located  on  the  axis,  the  solid  angle,  D,  can  be  calculated 
using  Equation  3. 


(3) 


In  Equation  3,  d  is  the  distance  front  the  source  to  the  detector  and  a  is  the 
radius  of  the  detector.  The  other  component  of  absolute  efficiency  is  the  detectors 
intrinsic  efficiency.  Intrinsic  efficiency  takes  into  account  how  many  of  the  neutrons 
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that  make  it  to  the  detector  are  actually  detected  by  the  detector  in  question.  It  can 
be  calculated  using  Equation  4, 


tint 


NHCTH 


Nhcth  +  Ncctc 


^  _  g—(NHo'H+Nccrc)d'^ 


(4) 


In  Equation  4,  N  is  the  number  density  of  atoms  in  the  detector  per  cm3,  a  is  the 
scattering  cross  section  of  the  given  element,  and  d  is  the  distance  traveled  by  the 
neutron  [24],  Once  the  intrinsic  efficiency  and  the  geometry  are  known,  the  absolute 
efficiency  can  be  calculated  using  Equation  5. 
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2.6  Muon  Detection 

Muon  detection  is  a  vital  component  to  any  muon  imaging  system  and,  therefore, 
must  be  considered  in  depth.  To  properly  associate  a  given  neutron  detection  with  an 
incident  muon,  the  muon  must  first  be  detected  and  verified  as  a  muon.  It  has  been 
proposed  to  utilize  a  series  of  Nal(Tl)  scintillation  detectors  to  accomplish  this  task  [1] . 
In  such  a  method,  the  incident  particle  will  deposit  energy  in  the  detectors  which  can 
then  be  utilized  to  both  confirm  the  presence  of  an  incoming  muon  and  provide  a 
means  of  energy  determination.  Muons  with  energies  greater  than  several  keV  will 
lose  a  portion  of  their  energy  within  the  detectors  through  ionization  processes  and 
produce  a  linear  energy  transfer  characterized  by  the  Bethe-Block  equation  (Equation 
6)  [25], 
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a;  =  log10  — (9) 

TTlpartC 

The  variable  re  is  the  electron  radius,  2.187  x  1013  m;  me  is  the  electron  mass, 
9.109  x  10~31  kg;  mpart  is  the  particle  mass;  Na  is  Avagadro’s  number  of  6.022  x  1023 
atoms/mol;  /  is  the  mean  ionization  potential  in  MeV;  Z  is  the  proton  number;  A 
is  the  atomic  mass;  p  is  the  material  density;  f3  =  v/c  where  v  is  the  velocity  of  the 
particle  and  c  is  the  speed  of  light;  7  is  the  relativistic  constant;  and  p  is  the  particle’s 
momentum.  In  each  equation  the  energy  is  supplied  in  MeV  and  masses  are  supplied 
in  MeV/c2.  C,a,x  1,  and  k  are  all  unitless  fitting  parameters  for  a  given  material 
which  have  been  tabulated  [26].  The  final  curve  for  a  muon  incident  upon  a  piece  of 
copper  material  is  shown  in  Figure  8. 

The  Bethe-Block  equation  can  be  used  in  conjunction  with  Equation  10,  where 
p  is  the  material  density  and  R  is  the  distance  traveled  in  the  material,  to  find  the 
energy  deposited  within  a  certain  material.  These  calculations  can  be  applied  to 
find  an  incident  muon  energy  by  using  a  detector  that  has  the  capability  to  measure 
deposited  energy  within  the  material. 
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Such  a  case  applies  for  utilizing  Nal(Tl)  detectors.  A  muon’s  energy  deposition 
curve  has  been  plotted  in  Figure  9  for  negative  muon  in  a  Nal(Tl)  detector.  As  can 
be  seen,  the  energy  deposition  varies  with  incident  muon  energy,  and  thus,  the  muon 
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Figure  8.  The  plot  displays  the  linear  energy  transfer  for  muons  passing  through  copper 
material  for  various  incident  muon  energies.  Copyright  (2001)  Academic  Press  [25] 


energy  can  be  calculated  based  upon  the  detector  response.  By  only  considering  events 
that  deposit  energies  of  anticipated  values,  many  sources  of  background  radiation  may 
be  eliminated.  This  method  can  then  be  utilized  to  verify  that  the  detectors  have 
detected  a  muon  instead  of  other  sources. 
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Incident  Muon  Energy  [MeV] 

Figure  9.  The  plot  displays  the  energy  deposited  in  a  Nal(Tl)  scintillation  detector 
based  on  the  incident  muon  energy.  All  data  has  been  calculated  according  to  the 
Bethe-Block  equation.  Reproduced  from  Dr.  Van  Dyk’s  dissertation.  [1] 
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III.  Methodology 


3.1  Overview 

In  order  to  conduct  this  experiment,  a  setup  was  created  in  which  incoming  muons 
were  detected  in  coincidence  with  neutrons  produced  by  muon  interactions  within  two 
materials,  iron  and  lead.  To  accomplish  such  a  task,  an  experimental  design  laid  out 
by  Dr.  Van  Dyk  was  utilized.  The  proposed  setup  featured  a  novel  muon  funnel  that 
directs  cosmic  muons  through  various  scattering  angles  and  has  been  shown  to  in¬ 
crease  the  muon  flux  by  3%  through  a  given  area  [1].  Immediately  below  this  funnel, 
a  series  of  four  Nal(Tl)  scintillation  detectors  were  used  to  determine  the  presence  of 
a  muon  and  it’s  incident  energy  on  the  system.  After  passing  through  the  Nal(Tl) 
detectors,  the  muon  then  entered  the  desired  material  to  be  studied.  The  experiment 
was  surrounded  with  24  EJ-309  pulse  shape  discriminating  liquid  scintillators.  These 
detectors  provided  a  method  to  detect  the  desired  neutron  products  and  also  allowed 
pulse  shape  discrimination  to  be  applied  to  the  signal  in  order  to  ignore  any  extrane¬ 
ous  gamma  signals  that  were  detected.  The  entire  detection  setup  below  the  funnel 
assembly  can  be  seen  in  Figure  10. 

3.2  Experimental  Setup 

Muon  Detection. 

The  muon  detection  system  served  to  both  detect  incoming  muons  in  coincidence 
with  the  produced  neutrons,  and  also  to  determine  the  incident  muon  energy  to 
analyze  any  energy  dependence  of  the  neutron  output.  To  accomplish  this  objective, 
four  Nal(Tl)  scintillation  detectors,  5.08  cm  thick,  40.64  cm  long,  and  10.16  cm  wide, 
were  located  directly  beneath  the  muon  funnel.  They  were  arranged  in  to  provide 
maximum  probability  of  a  vertical  muon  path  with  only  a  10°  entrance  window.  Each 
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Figure  10.  This  illustration  shows  the  A)  Nal(Tl)  scintillation  detectors,  B)  logic 
coincidence  trigger,  and  the  C)  neutron  detector  array  used  to  detect  cosmic  muons 
in  coincidence  with  their  produced  neutrons  within  the  material  to  be  examined.  The 
material  was  placed  at  the  center  of  the  neutron  detector  array. 


Nal(Tl)  detector  was  powered  using  an  Ortec  556  high  voltage  power  supply  at  1200 
Volts.  The  output  of  each  detector  went  directly  into  an  Ortec  113  preamplifier  with 
a  200  picofarad  capacitance  selected.  The  detector-preamplifier  combinations  were 
used  primarily  to  determine  a  coincident  muon  event,  and,  as  a  secondary  objective, 
to  determine  the  energy  spectrum  of  the  incident  muons.  To  accomplish  both  tasks 
simultaneously,  the  preamplifier  output  was  split.  The  primary  signal  was  passed  on 
to  an  Ortec  935  quad  constant  fraction  discriminator  and  then  to  an  Ortec  4020  quad 
4-input  logic  unit.  The  logic  unit  was  set  to  output  a  signal  only  in  the  case  where  all 
four  Nal(Tl)  received  a  pulse  simultaneously.  In  such  a  case,  the  output  was  used  as 
the  coincident  trigger  and  was  fed  to  the  external  trigger  input  for  the  CAEN  V1724 
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and  three  CAEN  V1720  digitizers  with  a  CAEN  V1718  controller.  The  secondary 
signal  was  fed  directly  into  the  input  channels  0,  1,  2,  and  3  of  the  CAEN  1724 
board.  This  data  was  filtered  and  captured  by  the  DPP-PHA  pulse  height  analysis 
software  provided  by  CAEN  with  settings  seen  in  Appendix  A.  The  output  was  saved 
using  list  mode  and  post  processed  through  a  Matlab  script  which  can  be  seen  in 
Appendix  E.  The  muon  detection  portion  of  the  experimental  setup  is  illustrated  in 
Figure  11. 


Figure  11.  This  illustration  depicts  the  A)  muon  funnel,  B)  Nal(Tl)  scintillation  de¬ 
tectors,  and  C)  logic  coincidence  trigger  used  to  determine  the  initial  muon  entrance 
into  the  material. 
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Neutron  Detection. 


Beneath  the  Nal(Tl)  scintillation  detectors,  an  array  of  24  EJ-309  liquid  scintilla¬ 
tion  detectors  with  five  inch  photomultiplier  tubes  were  arranged  into  four  rings  of  six 
detectors  apiece.  One  such  ring  can  be  seen  in  Figure  12,  while  the  entire  array  can 
be  seen  in  Figure  13.  All  the  detectors  were  powered  through  a  CAEN  SY4527  crate 
with  a  24  channel  3  kV  power  supply  module.  Each  detector  had  an  optimum  applied 
voltage  level  given  by  the  manufacturer,  Eljen  Technologies.  These  voltages  are  listed 
in  Appendix  C  for  reference.  The  detector  outputs  were  passed  to  CAEN  V1720 
digitizers  where  the  signals  were  recorded  and  PSD  values  were  calculated  using  the 
DPP-PSD  firmware  provided  by  CAEN.  The  settings  for  the  DPP-PSD  firmware  can 
be  seen  in  Appendix  B. 


Figure  12.  The  EJ-309  neutron  detectors  were  constructed  to  fit  six  per  ring  of  the 
neutron  array  with  enough  room  in  the  center  to  place  material  for  inspection. 

Previously  a  77  mCi  plutonium  beryllium  source  had  been  utilized  to  examine 
the  PSD  values  produced  by  the  liquid  scintillators  used  in  this  experiment.  Figure 


22 


Figure  13.  The  neutron  detector  array  is  composed  of  24  EJ-309  liquid  scintillation 
detectors.  This  array  provides  nearly  a  Ai r  solid  angle  of  detection  capability  when 
material  is  placed  directly  at  the  center  of  the  array  (A). 
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14  displays  the  acquired  2D  histogram  for  both  the  gammas  and  neutrons  that  were 
detected  by  the  liquid  scintillators.  A  filter  was  applied  to  the  PSD  values  in  post 
processing  to  eliminate  as  many  gamma  events  as  possible  by  choosing  cut  lines  at 
PSD  values  of  0.1  and  0.3  for  the  lower  and  upper  limits  respectively. 
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Figure  14.  PSD  2D  histogram  used  to  determine  neutron  discrimination  cut  lines. 
The  lower  and  upper  cut  lines  were  selected  at  PSD  values  of  0.1  and  0.3  respectively. 
Reproduced  from  Dr.  Van  Dyk’s  dissertation.  [1] 

Digitizer  Synchronization. 

An  incompatibility  in  the  clock  synchronization  between  the  two  different  CAEN 
digitizer  boards,  V1720  and  V1724,  prevented  accurate  correlation  between  the  inci¬ 
dent  muon  events  and  the  neutrons  that  were  produced.  To  solve  this  problem,  one 
channel  on  each  of  the  three  CAEN  VI 720  digitizers  was  switched  from  one  of  the 
liquid  scintillator  detectors  and  replaced  with  an  input  signal  from  a  Hewlett  Packard 
33120A  function  generator  which  can  be  seen  in  Figure  15.  A  1  V,  10  MHz  square 
pulse  was  continuously  provided  to  each  digitizer  to  ensure  a  recorded  time  stamp  ev¬ 
ery  time  the  external  trigger  was  set  off  by  an  incident  muon.  In  order  to  accomplish 
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this,  three  of  the  bottom  neutron  detectors  were  taken  offline  and  their  data  was  not 
gathered.  Once  the  data  had  been  gathered  for  each  run,  a  Matlab  script  was  used 
to  correlate  the  timestamps  for  each  neutron  event  to  one  incident  muon.  This  script 
can  be  seen  in  Appendix  E. 


Figure  15.  A  function  generator  was  used  to  ensure  time  stamps  for  every  neutron 
detection  were  recorded  with  each  muon  event,  since  clock  synchronization  between 
the  two  different  digitizer  boards  was  not  possible. 
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IV.  Results  and  Analysis 


4.1  Modeling  Results 

To  examine  the  feasibility  of  using  the  phenomenon  of  //IF  as  a  detection  method, 
the  muon  induced  neutron  production  rates  need  to  be  determined  for  muon  interac¬ 
tions  with  fissile  materials.  A  simulation  was  performed  to  examine  this  characteristic 
utilizing  the  Monte  Carlo  simulation  tool  Geant4.  This  tool  is  produced  by  CERN 
to  simulate  the  interaction  and  passage  of  high  energy  particles  through  matter  and 
can  be  easily  utilized  to  perform  the  necessary  analysis. 

Initial  Modeling. 

Initially,  a  5  cm  radius  and  20  cm  tall  cylinder  was  simulated  as  the  target  material. 
A  47 r  spherical  detector  was  created  around  the  target  material  to  detect  and  measure 
the  emitted  neutron  energy  spectra  and  counts.  The  entire  volume,  excluding  the 
target  material,  was  filled  with  air  composed  of  70%  nitrogen  and  30%  oxygen.  Muons 
were  simulated  as  an  isotropic  point  source  located  at  the  origin  and  centered  within 
the  target  material  allowing  for  maximum  interaction.  The  entire  geometry  setup  can 
be  seen  in  Figure  16. 

Five  different  target  materials  were  used  representing  a  wide  range  of  atomic 
numbers  and  enrichment  levels.  These  materials  included  enriched  uranium,  with 
an  enrichment  level  of  90%235U  and  10%238U,  pure  238U,  lead,  iron,  and  aluminum. 
Each  of  these  materials  was  exposed  to  100,000  muons  at  one  of  two  different  initial 
energies.  The  first  trial  was  conducted  at  a  muon  source  energy  of  1  GeV.  This  energy 
was  selected  as  it  is  on  the  order  of  the  average  energy  of  the  cosmic  muon  flux  at 
ground  level  after  being  moderated  by  30  cm  of  concrete  as  seen  in  Figure  17.  A 
second  trial  was  conducted  with  a  muon  source  energy  of  100  MeV.  This  energy  was 


26 


ma 


Figure  16.  The  geometry  setup  for  the  Geant4  simulation  modeling  muon  induced 
neutrons.  This  figure  does  not  display  the  4n  detector  surrounding  the  target  material. 
The  green,  red,  and  blue  tracks  represent  neutral,  negative,  and  positively  charged 
particles  created  from  the  initial  muon  interactions. 

selected  as  it  is  on  the  lower  side  of  the  cosmic  muon  spectrum,  but  has  much  higher 
neutron  production  rate  within  the  target  materials. 

100  MeV  Muon  Source. 

As  can  be  seen  from  Figure  18,  the  neutron  energy  spectra  from  each  of  four 
materials  have  both  quantitative  and  qualitative  differences.  As  expected  the  enriched 
uranium  target  produced  the  most  neutrons  with  the  238U  target  following  closely 
behind.  However,  both  uranium  targets  produced  50%  more  neutrons  than  the  lead 
target  and  a  factor  of  3  greater  than  the  iron  target.  It  should  also  be  mentioned 
that  although  a  trial  run  was  conducted  for  each  of  the  five  materials,  aluminum  did 
not  produce  a  statistically  significant  number  of  neutrons  to  be  able  to  construct  a 
spectrum. 

1  GeV  Muon  Source. 

As  the  energy  of  the  muon  source  was  increased  to  1  GeV,  the  number  of  neutrons 
produced  were  reduced.  In  fact,  although  the  increase  in  energy  was  only  one  order 
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Figure  17.  This  plot  displays  the  muon  energy  spectrum  after  bing  moderated  by  30 
cm  of  concrete,  much  like  has  been  done  in  this  experiment.  Secondary  neutrons  are 
created  through  muon  interactions  with  the  concrete.  The  figure  has  been  reproduced 
from  work  in  the  public  domain  by  the  Department  of  Energy.  [27] 
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Figure  18.  Muon  induced  fission  neutron  energy  spectrums  from  four  different  materials 
produced  using  100  MeV  muons.  It  should  be  noted  that  although  aluminum  was  also 
simulated,  the  neutron  production  was  negligible  and  no  spectrum  was  produced. 


of  magnitude,  the  neutron  count  decreased  two  orders  of  magnitude  in  uranium  and 
three  orders  of  magnitude  in  lead  and  iron.  In  addition  to  the  lack  of  counts,  the 
spectra  became  much  more  sparse  than  their  counterparts  corresponding  to  the  same 
material  with  the  100  MeV  source  but  reducing  or  removing  any  noticeable  features. 
These  spectrums  can  be  seen  in  Figure  19.  Once  again,  the  simulation  ran  all  five 
materials.  However,  similar  to  the  100  MeV  muon  source,  the  aluminum  target  did 
not  produce  a  single  neutron  and  is  therefore  not  displayed. 
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Figure  19.  Muon  induced  fission  neutron  energy  spectrums  from  four  different  materials 
produced  using  1  GeV  muons.  It  should  be  noted  that  although  aluminum  was  also 
simulated,  the  neutron  production  was  negligible  and  no  spectrum  was  produced. 


Neutron  Yield  per  Muon. 

When  considering  the  detection  capability  and  usefulness  of  muon  induced  fission, 
it  is  important  to  keep  in  mind  that  the  neutrons  produced  from  each  fission  event 
may  act  as  a  secondary  signal  to  be  acquired  and  analyzed.  Each  additional  signal 
gathered  has  the  potential  to  deliver  a  greater  confidence  of  the  presence  of  SNM. 
Table  1  shows  the  average  number  of  neutrons  produced  for  each  muon  simulated.  As 
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expected,  the  enriched  uranium  produced  the  most  neutrons  while  aluminum  yielded 
the  least. 


Table  1.  Average  Neuton  Yield  per  Muon  For  Initial  Simulation 


100 

MeV 

1 

GeV 

Enriched  Uranium 

2.3 

± 

0.1 

0.1 

± 

0.1 

Depleted  Uranium 

1.5 

± 

0.1 

0.1 

± 

0.1 

Lead 

1.3 

± 

0.1 

0.1 

± 

0.1 

Iron 

0.8 

± 

0.1 

0.1 

± 

0.1 

Final  Modeling. 

After  the  initial  modeling  attempt,  a  more  accurate  simulation  of  the  experiment 
was  required.  For  this  simulation,  a  block  of  material  was  created  as  the  target  and 
located  at  the  center  of  the  world  volume.  A  47t  spherical  detector  with  a  1  meter 
radius  was  created  around  the  target  material  to  detect  the  number  of  neutrons  pro¬ 
duced.  The  entire  volume,  excluding  the  target  material,  was  filled  with  air  composed 
of  70%  nitrogen  and  30%  oxygen.  Muons  were  created  at  a  distance  of  1.5  meters  from 
the  origin  and  given  an  initial  velocity  towards  the  target  block.  Upon  interaction, 
only  immediate  daughter  products  of  muon  interactions  were  found  and  recorded,  as 
well  as  any  neutrons  that  passed  through  the  detection  sphere.  In  these  scenarios, 
three  materials  were  selected  to  be  modeled,  90%  enriched  uranium,  lead,  and  iron. 
Each  material  had  a  20  cm  x  10  cm  surface  exposed  to  the  incident  muons,  with 
the  thicknesses  varied.  The  lead  and  enriched  uranium  were  5  cm  thick,  while  the 
iron  and  another  lead  simulation  were  conducted  with  a  target  block  thickness  of  15 
cm.  These  thicknesses  were  selected  as  they  corresponded  to  the  available  material 
on  hand  for  the  physical  experiment.  An  example  of  the  target  material  and  its  setup 
can  be  seen  in  Figure  20. 
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Figure  20.  The  geometry  setup  for  a  Geant4  simulation  modeling  muon  induced  neu¬ 
trons  in  a  5  cm  thick  lead  target.  This  figure  does  not  display  the  47r  detector  surround¬ 
ing  the  target  material.  The  green  tracks  represent  neutral  particles  created  from  the 
initial  muon  interactions,  while  the  red  track  is  the  incident  muon  path. 

To  gain  a  better  understanding  of  the  muon  interactions  at  higher  energies,  two 
additional  energy  simulations  were  conducted.  Initial  muon  energies  of  100  MeV,  1 
GeV,  10  GeV,  and  100  GeV  were  selected  and  100,000  muons  were  simulated  for  each 
energy  on  each  material.  The  resulting  neutron  yields  and  their  errors  can  be  seen 
plotted  in  Figure  21. 

As  can  be  seen  in  Figure  21,  the  100  MeV  incident  muons  produced  several  orders 
of  magnitude  more  average  number  of  neutrons  than  those  at  higher  energies.  How¬ 
ever,  there  is  an  increase  in  neutron  production  with  an  increase  in  incident  muon 
energy  greater  than  1  GeV.  Inspection  into  the  documentation  on  the  cross  sections 
utilized  by  Geant4  revealed  that  little  is  known  about  the  cross  sections  in  the  region 
of  low  muon  energy  below  approximately  1  GeV  for  neutron  production.  Inelastic 
interactions  between  muons  and  nuclei  gains  importance  at  energies  above  10  GeV, 
while  below  this  threshold  minimal  neutron  production  occurs  via  this  method.  How- 
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Figure  21.  The  dependence  of  average  neutron  yield  on  muon  energy  for  four  different 
material  and  geometry  configurations. 

ever,  in  the  event  that  the  muon  can  be  thermalized,  the  capture  process  dominates 
within  the  simulations  though  the  Geant4  process  name  of  muMinusCaptureAtRest, 
but  no  cross  sections  are  listed.  This  is  indeed  in  line  with  the  neutron  yields  pro¬ 
duced  within  this  simulation.  Total  neutron  yields  of  the  final  simulation  can  be  seen 
in  Table  2. 


Table  2.  Average  Neuton  Yield  per  Muon  For  Final  Simulation 


100  MeV 

1  GeV 

10  GeV 

100  GeV 

5  cm  Uranium 

6.6  ±  0.1 

0.1  ±  0.1 

0.1  ±  0.1 

0.1  ±  0.1 

5  cm  Lead 

1.4  ±  0.1 

0.1  ±  0.1 

0.1  ±  0.1 

0.1  ±  0.1 

15  cm  Lead 

1.4  ±  0.1 

0.1  ±  0.1 

0.1  ±  0.1 

0.1  ±  0.1 

15  cm  Iron 

1.3  ±  0.1 

0.1  ±  0.1 

0.1  ±  0.1 

0.1  ±  0.1 

In  terms  of  detection  capabilities,  the  neutron  yield  at  the  lower  muon  energies  may 
be  able  to  be  utilized  as  a  secondary  detection  method.  However,  when  considering 
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atmospheric  muon  energies,  simulations  suggest  most  muons  will  not  contribute  to 
the  generation  of  this  signal  unless  they  are  previously  moderated  to  lower  energies. 
If  it  could  be  arranged  for  some  type  of  moderating  material  to  be  placed  between 
the  target  and  the  source,  an  increase  in  neutron  production  may  be  seen. 

4.2  Experimental  Results 

To  perform  the  experiment,  blocks  of  iron  and  lead  were  placed  at  the  center  of  the 
array  of  neutron  detectors.  These  materials  were  chosen  because  of  their  availability 
and  varied  atomic  number.  Each  material  had  an  exposed  area  of  20  cm  x  10  cm 
exposed  to  the  incident  muons,  but  the  thickness  was  varied.  The  iron  block  had  a 
thickness  of  15  cm  while  two  blocks  of  lead  were  used  with  thicknesses  of  5  cm  and  15 
cm.  Every  trial  was  performed  for  a  duration  of  seven  days  with  an  additional  seven 
day  measurement  of  natural  background  neutron  levels. 

Muon  Energy  Determination. 

In  order  to  properly  classify  an  event  as  muon-induced,  the  muon  event  was  de¬ 
tected  in  coincidence  by  four  Nal(Tl)  scintillation  detectors.  These  detectors  captured 
pulse  height  information  which  was  converted  to  energy  through  a  calibration  method 
which  is  detailed  in  Appendix  D.  The  Nal(Tl)  detectors  were  calibrated  using  known 
gamma  sources  of  Cs-137  and  Co-60.  In  addition  to  the  low  energy  gamma  sources, 
a  higher  energy  point  was  used  with  the  most  likely  energy  deposited  in  the  detec¬ 
tors  corresponding  to  the  minimum  ionizing  potential  for  the  Nal(Tl)  detectors  found 
using  Equation  6  to  be  24.297  MeV. 

Once  this  calibration  was  applied,  pulse  heights  could  be  converted  to  energy 
deposited  within  the  Nal(Tl)  detectors.  The  energy  deposited  was  then  filtered  to 
remove  any  energy  deposition  events  not  consistent  with  that  anticipated  by  the 
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Bethe-Block  equation  for  a  muon  event.  The  minimum  threshold  was  chosen  to  be 
the  minimum  ionizing  potential  of  24.297  MeV,  while  the  maximum  threshold  was 
chosen  to  be  the  expected  energy  deposition  of  a  5  GeV  muon  of  31.75  MeV.  Each 
threshold  was  also  modified  to  compensate  for  a  detector  resolution  of  12%  by  setting 
the  threshold  value  at  3 a  below  and  above  the  stated  values,  respectively. 

After  threshold  values  had  been  applied  for  each  muon  event,  a  muon  energy  was 
varied  from  140  -  5000  MeV  in  1  MeV  increments  for  a  muon  event  incident  upon 
the  first  Nal(Tl)  detector.  Each  of  these  incident  energy  values  was  then  used  to 
determine  the  theoretical  energy  deposition  in  each  of  the  four  Nal(Tl)  detectors 
using  Equation  6.  The  theoretical  values  were  then  compared  to  the  measured  values 
using  Pearson’s  y2  test  statistic  calculated  using  Equation  11. 


2  V''  ( dEiBB 

x  =2^ — — 

2=1 


d-Ey measured) 

dEiBB 


where  dEhBB  is  the  theoretically  predicted  energy  deposited  in  each  of  the  Nal(Tl) 
detectors  and  dEi  measured  is  the  experimentally  measured  energy  deposited. 

The  initial  muon  energy  that  resulted  in  the  lowest  \2  test  statistic  was  the  most 


likely  energy  of  the  muon  incident  upon  the  system.  The  theoretical  energy  of  the 
muon  leaving  the  last  Nal(Tl)  detector  was  then  calculated  and  used  as  the  muon 
energy  incident  upon  the  target.  The  final  muon  energy  distribution  can  be  seen  in 


Figure  22.  Previous  experiments  utilizing  the  same  experimental  setup  have  found 
the  mode  and  mean  of  the  distribution  to  be  550  MeV  and  926  MeV  respectively  [1] . 
However,  this  series  of  experiments  found  the  mode  and  mean  of  the  distribution  to 


be  approximately  180  MeV  and  520  MeV  respectively.  The  discrepancy  likely  rises 
from  a  change  in  the  operating  voltage  of  the  Nal(Tl)  detectors  and  differences  in 
calibration  methods. 
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Figure  22.  This  figure  displays  the  incident  muon  energy  distribution  after  passing 
through  the  muon  funnel  and  all  four  Nal(Tl)  detectors.  The  slight  spike  located  at 
approximately  200  MeV  corresponds  to  the  location  of  the  minimum  ionizing  potential 
in  the  Bethe-Block  equation. 


Neutron  Yield  Analysis. 

The  number  of  neutrons  produced  for  each  muon  event  were  recorded  and  are 
displayed  in  Figure  23. 

As  was  expected,  when  no  material  was  present  the  majority  of  incident  muons 
were  detected  in  coincidence  with  zero  neutrons.  The  number  of  muons  producing  one 
neutron  decreased  by  an  order  of  magnitude,  and  again  by  another  order  of  magnitude 
when  two  neutrons  were  produced.  However,  an  interesting  feature  arose  when  the 
neutron  counts  were  observed  with  material  present.  In  all  but  one  case,  the  number 
of  muon  events  that  produced  zero  neutrons  was  less  than  the  amount  that  produced 
one  neutron.  This  suggests  that  a  majority  of  the  muons  that  are  incident  upon  the 
material,  interact  in  such  a  way  as  to  produce  at  least  one  neutron  and  is  supported  by 
the  average  neutron  calculation  based  upon  the  simulation  results.  This  assumption 
is  supported  further  by  examining  the  neutron  counts  once  the  background  neutron 
levels  have  been  removed  as  seen  in  Figure  24.  In  such  a  case,  a  change  in  the 
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Figure  23.  This  figure  displays  the  frequency  of  occurrence  for  each  neutron  yield  over 
the  course  of  the  seven  day  experiment. 

maximum  neutrons  produced  is  seen  as  an  increase  from  one  neutron  being  produced 
for  an  incident  muon  to  two  produced  neutrons. 

In  order  to  directly  compare  the  neutron  yields  from  the  different  materials  and 
geometries,  the  muon  counts  were  scaled  to  produce  the  probability  distribution  dis¬ 
played  in  Figure  25.  These  calculated  probability  distributions  display  the  likelihood 
that  a  specific  number  of  neutrons  will  be  detected  above  the  background  levels  based 
upon  the  data  gathered  over  the  seven  day  period  for  each  sample.  A  two  sample 
t-test  was  performed  comparing  each  data  set  to  the  other  two,  and  the  null  hypoth¬ 
esis  of  ji\  =  /i2  could  not  be  rejected  in  any  case  with  the  lowest  p-value  being  86% 
suggesting  that  all  data  sets  are  likely  correlated. 

Mean  neutron  yield  per  detected  muon  was  calculated  for  the  experimental  results 
in  the  same  manner  as  the  simulated  results,  using  total  muons  and  total  neutrons 
detected.  These  results  displayed  in  Table  3.  Immediately  an  anomaly  arose  when 
inspecting  the  average  neutron  yield  per  muon,  as  the  15  cm  block  of  lead  produced 
the  lowest  average  neutron  yield.  However,  for  a  period  of  time  the  power  to  the 
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Neutron  Counts  per  Muon 

Figure  24.  This  figure  displays  the  frequency  of  occurrence  for  each  neutron  yield  after 
background  levels  were  removed. 

neutron  detectors  was  temporarily  shut  off  while  power  to  the  muon  coincidence  unit 
remained  functional.  The  result  of  this  power  fluctuation  resulted  in  a  higher  number 
of  muons  passing  through  the  material  than  the  neutron  data  accounts  for.  This 
difference  also  explains  the  uniqueness  of  the  15  cm  lead  data  in  Figure  23  having  a 
higher  muon  count  at  a  neutron  yield  of  zero  than  of  one. 


Table  3.  Mean  neutron  yields  per  muon  after  7  days  of  exposure 


Material 

Total 

Above  Background 

100  MeV 

1  GeV 

5  cm  Lead 

2.8  ±  0.1 

4.0  ±  0.1 

2.1  ±  0.1 

2.1  ±  0.3 

15  cm  Lead 

2.2  ±  0.1 

3.8  ±  0.1 

1.6  ±  0.1 

1.1  ±  0.2 

15  cm  Iron 

3.4  ±  0.1 

4.8  ±  0.1 

3.0  ±  0.2 

2.1  ±  0.3 

Additionally,  the  mean  values  for  the  number  of  neutrons  produced  above  the 
background  levels  were  calculated  and  can  also  be  found  in  Table  3. 


37 


0.25 


— 15  cm  Iron 


Figure  25.  The  probability  distribution  for  an  incident  muon  to  produce  a  given  neutron 
count  above  background  levels. 

Energy  Dependent  Neutron  Yield. 

To  further  compare  the  neutron  output  found  experimentally,  the  neutron  yields 
were  plotted  against  incident  muon  energy,  and  can  be  seen  in  Figure  26.  In  these 
plots,  each  point  represents  one  muon  event.  A  qualitative  increase  in  neutron  pro¬ 
duction  can  be  seen  visually  between  the  background  measurement  and  the  addition 
of  any  material. 

To  further  analyze  this  information,  the  muon  events  were  binned  by  energy  in  100 
MeV  increments  and  their  average  neutron  yields  calculated.  At  lower  energies  a  more 
consistent  neutron  yield  was  produced,  but  the  noise  increased  with  incident  muon 
energy.  To  gain  a  better  understanding  of  this  phenomena,  the  neutron  yield  was 
plotted  with  the  corresponding  error  in  each  point  in  Figure  27.  The  increased  noise 
at  higher  incident  muon  energies  is  to  be  expected  as  the  muon  energy  distribution 
was  primarily  below  1000  MeV  leading  to  very  poor  counting  statistics  in  the  higher 
energy  region. 

Another  two  sample  t-test  was  conducted  to  compare  each  of  the  materials  average 
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Figure  26.  This  figure  displays  the  number  of  neutrons  produced  by  and  energy  of  each 
muon  incident  on  the  system.  Visual  inspection  reveals  that  neutron  counts  are  above 
background  level,  but  more  quantitative  analysis  must  be  done  to  distinguish  between 
materials. 


neutron  output  to  the  other  three  sets  of  data.  In  each  case  the  null  hypothesis  of 
hi  =  h 2  could  be  rejected  with  a  95%  probability.  This  result  suggests  that  the 
muon  energy  dependent  neutron  yields  of  each  material  were  not  equivalent.  It  is 
predicted  that  a  poor  PSD  filter  has  been  applied  to  the  data.  This  would  mean  that 
gamma  events  created  within  the  target  material  are  being  counted  as  neutrons  when 
detected  by  the  liquid  scintillation  detectors,  artificially  increasing  counts.  In  such  a 
case,  additional  gamma  events  would  reach  the  detector  when  produced  in  the  iron 
target  that  would  be  attenuated  when  created  within  the  lead  target. 

Finally,  the  incident  muon  energy  spectra  were  examined  based  upon  the  neutron 
yield  of  that  particular  event.  As  seen  in  Figure  28,  the  spectra  shape  are  similar  to 
the  incident  muon  energy  spectrum  seen  in  Figure  22.  Also  important  to  note  are  the 
relative  heights  of  each  material’s  spectra.  With  a  neutron  yield  of  zero,  one,  or  two 
neutrons  created,  the  order  is  as  expected.  However,  when  three  or  more  neutrons 
are  produced,  the  peaks  of  the  spectra  begin  to  overlap  and  become  indistinguishable 
from  one  another.  These  features  agree  with  the  previous  findings  seen  in  Figure 
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Figure  27.  This  plot  displays  the  average  neutrons  produced  for  a  muon  with  a  given 
energy  incident  upon  the  system. 


25  with  both  lead  samples  yielding  more  events  with  a  neutron  yield  of  one  and  two 
neutrons,  but  the  overlap  between  lead  and  iron  begins  to  occur  with  a  three  detected 
neutrons. 
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Figure  28.  This  figure  displays  the  muon  energy  spectra  for  various  neutron  yields. 
As  expected  it  follows  the  shape  of  the  incident  muon  energy  distribution  in  all  cases. 
Above  a  neutron  yield  of  three  neutrons  the  muon  counts  continue  to  decrease,  but  the 
general  shape  of  the  spectrum  remains  the  same. 


41 


V.  Conclusions 


5.1  Research  Summary 

As  the  need  for  better  imaging  and  noninvasive  interrogation  techniques  rises, 
the  applicability  of  muon  imaging  techniques  grows.  In  these  scenarios,  cosmic  and 
atmospheric  muons  can  offer  superior  penetrating  capabilities  and  can  be  utilized 
to  determine  the  contents  of  a  sealed  container  such  as  a  shipping  container  or  a 
nuclear  weapon.  However,  due  to  the  requirement  of  a  detector  on  both  sides  of  the 
target  material,  the  feasibility  of  utilizing  these  techniques  for  real  world  applications 
diminishes.  Secondary  signals  may  help  to  provide  an  alternate  means  of  utilizing 
muon  imaging  in  telescopic  mode.  Initial  findings  suggest  neutrons  produced  by 
muon  induced  fission  may  be  able  to  provide  one  such  secondary  signal. 

Initial  simulations  of  two  different  muon  energies,  100  MeV  and  1  GeV,  were 
performed  on  five  different  materials  of  various  atomic  numbers.  Neutron  energy 
spectrums  were  gathered  and  analyzed  for  each  of  the  five  materials.  The  neutron 
yields  per  muon  were  calculated  and  found  to  range  from  between  2.3  ±  0.1  for 
enriched  uranium  down  to  negligible  amounts  for  aluminum  when  exposed  to  the 
100  MeV  source.  As  the  energy  of  the  muons  was  increased  to  1  GeV,  the  neutron 
yields  shrank  to  negligible  levels.  These  simulations  also  suggested  that  there  is 
little  difference  in  neutron  yield  produced  in  non-fissile  material  while  the  neutrons 
produced  in  an  enriched  uranium  block  was  greater  by  a  factor  of  five. 

Experimental  results  differed  from  these  findings.  A  probability  distribution  was 
constructed  from  the  neutron  yields  of  each  incident  muon.  The  various  probability 
distributions  produced  by  both  the  iron  and  lead,  as  well  as  the  two  different  thick¬ 
nesses  of  lead  were  nearly  identical.  Statistical  T-tests  were  conducted  to  compare 
the  various  distributions  to  one  another.  In  no  case  were  the  results  conclusive  and 
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able  to  reject  the  null  of  ji\  =  /i2  suggesting  minimal  difference  in  the  distributions. 
The  average  neutron  yield  was  also  calculated  and  found  to  be  3.4  ±  0.1  for  a  15  cm 
thick  block  of  iron,  2.8  ±  0.1  for  a  5  cm  thick  block  of  lead,  and  2.2  ±  0.1  for  a  15 
cm  thick  block  of  lead. 

The  energy  distribution  was  found  for  a  muon  incident  upon  the  system  and  av¬ 
erage  neutron  yields  were  calculated  for  each  100  MeV  region.  The  average  neutron 
yields  were  plotted  against  incident  muon  energy  and  statistical  T-tests  were  con¬ 
ducted  comparing  the  various  materials’  distributions  to  one  another.  In  every  case, 
the  null  hypothesis  of  /q  =  could  be  rejected  with  a  95%  probability  suggesting 
that  the  distributions  were  statistically  different.  Additionally,  incident  muon  en¬ 
ergy  spectra  were  examined  for  each  neutron  yield  output  and  differences  began  to 
diminish  above  three  detected  neutrons. 

From  these  findings,  muon  imaging  operated  in  telescopic  mode  may  be  able  to 
detect  SNM  within  an  enclosed  container  by  observing  the  neutron  output.  Simu¬ 
lations  indicate  that  highly  enriched  uranium  will  provide  a  significant  increase  in 
neutron  output  when  compared  to  other  materials  such  as  lead  and  iron.  In  addi¬ 
tion,  if  incident  muon  energy  can  be  determined  other  materials  may  be  able  to  be 
distinguished  by  examining  the  muon  energy  dependent  neutron  output. 

5.2  Future  Work 

Since  muon  imaging  holds  such  promise  for  accurate  detection  of  SNM  in  both 
portal  inspection  points  and  international  treaty  verification,  further  work  in  this  area 
is  highly  recommended.  As  a  follow  on  to  this  experiment,  a  more  accurate  method  of 
muon  energy  determination  would  benefit  the  accuracy  of  any  findings.  Another  topic 
of  study  which  would  be  helpful  in  the  determination  of  neutron  output,  would  be  to 
examine  in  further  detail  the  PSD  values  for  neutron  and  gamma  events  in  each  liquid 
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scintillator  as  well  as  the  PSD  value  of  a  muon  within  the  same  detectors  used  in  this 


experiment.  It  is  possible  that  atmospheric  muons  may  be  able  to  impart  a  similar 
PSD  value  as  the  neutrons  increasing  their  count  rate  and  thus  yield.  Finally,  the  most 
important  experiment  that  must  be  performed  is  the  analysis  of  muon  induced  neutron 
output  in  fissionable  elements.  While  simulations  do  suggest  neutron  production  at 
elevated  levels,  experimental  confirmation  is  needed. 
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Appendix  A.  PHA  Software  Settings 


Figure  29.  Computer  software  settings  for  pulse  height  analysis  of  the  Nal(Tl)  detec¬ 
tors. 
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Appendix  B.  PSD  Software  Settings 


Figure  30.  Computer  software  settings  for  pulse  shape  discrimination  and  neutron 
detection. 
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Appendix  C.  Neutron  Detector  Voltages 


Table  4.  Voltage  settings  for  the  various  liquid  scintillation  neutron  detectors. 


LVtrctor  Number 

Smal  Number 

lUgfc.  Noftjgr  ! 

1.1 

4567-01-21 
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1.2 
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1650 

1  J 
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i'.'iT-ii  i-^n 
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IS 
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:  r. 
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2.1 

4567*01-05 

1300 

~T2 

. 

1300 

1  3 

4>7  mil 

1250 

14 

02 

1300 

S) 

4567-01-06 

1400 

2  6 

4567-01-24 

1375 

3.1 

4567-01-04 
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i>7.t!|  a 

. 

US 

i>7  111-17 

130,1 

3-4 
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3.5 

;>7-Ul  16 
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1ft 

l>7  IJ1  [■ 
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l  1 

iii-ips 

1325 

1 
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.  ;j. 
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45 
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4567-01413 
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Appendix  D.  Nal(Tl)  Detector  Calibration 


To  properly  caliberate  the  Nal(Tl)  scintillation  detectors  two  gamma  sources  were 
selected  to  provide  known  energy  depositions  within  the  detectors.  A  Cs-137  and  a 
Co-60  source  were  selected  from  the  available  resources  and  provided  three  well  known 
gamma  energies  of  662  keV,  1173  keV,  and  1333  keV.  The  two  sources  were  exposed 
to  the  Nal(Tl)  detectors  and  spectra  were  gathered.  Each  of  these  spectra  were 
then  analyzed  and  the  observed  peak  values  were  recorded.  Because  muon  energy 
deposition  primarily  occurs  at  much  higher  energies  than  can  be  deposited  by  gamma 
sources,  an  additional  point  of  higher  energy  must  be  used  minimize  extrapolation 
error. 

To  more  accurately  predict  the  energy  deposition  at  these  higher  energies,  the 
minimum  ionizing  potential  can  utilized.  The  most  probable  muon  energy  deposi¬ 
tion  corresponds  to  the  minimum  ionizing  potential  of  22.297  MeV  in  the  Nal(Tl) 
detectors  according  to  the  Bethe-Block  equation,  Equation  6.  Upon  completion  of 
an  experimental  run,  the  pulse  height  spectra  were  plotted  for  each  of  the  Nal(Tl) 
detectors  and  can  be  seen  in  Figure  31.  The  maximum  value  of  these  peaks  were 
determined  to  represent  the  minimum  ionizing  potential  energy  as  has  been  done  in 
previous  experiments  [1]. 

Once  the  pulse  height  channels  were  gathered  from  the  minimum  ionizing  potential 
peak  and  the  known  gamma  sources,  they  were  plotted  against  the  known  energies. 
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Figure  31.  The  peak  value  of  the  middle  peak  in  the  pulse  height  distributions  were  de¬ 
termined  to  be  the  minimum  ionizing  potential  according  to  the  Bethe-Block  equation 
for  Nal(Tl)  with  an  energy  of  24.297  MeV. 


The  calibration  equations  were  calculated  to  be 

Detectorl  :  y  =  400x  +  220 
Detector2  :  y  =  380x  +  170 
Detector?*  :  y  =  360x  +  100 
Detector 4  :  y  =  450x  +  670 


and  can  be  seen  in  Figure  32. 
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Figure  32.  Nal(Tl)  calibration  curves  for  each  of  the  four  Nal(Tl)  detectors.  The 
calibration  points  with  the  associated  error  can  be  seen  as  well  at  0.662  MeV,  1.173 
MeV,  1.333  MeV,  and  24.297  MeV. 
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Appendix  E.  Matlab  Code 


9'2'2'2'2'2'9'2'9-9'9'2-2'2-2'2'2'2'9'9'9'9'2-2'2'2'2'9'9'2'9'2'9'2-2'2'2'2'2'9-9'9'2'2'2- 

ooooooooooooooooooooooooooooooooooooooooooooo 


9'2'9-2'9-2'9-9'9-9'9'9-2'9-2'9-2'9-9'9'9'2'9-2'9'2'9'9'9-9'2'9'9'9-9'2'2'9-9'9-9'2'9'9'9- 

ooooooooooooooooooooooooooooooooooooooooooooo 


% IMPORT  Data 


9'2'9'2'9-2'9-9'9-9'9'9-2'9-2'9-2'9-9'9'9'2'9-2'9'2'9-9'2-9'9'9'9'9-2'2'2'9-9'9-9'9'9'2'9- 

ooooooooooooooooooooooooooooooooooooooooooooo 


9'2'2-2'2'2'9-2'2-9'9'2-2'2-2'2'2'2'9'9'9'9'2-2'2-2'2'9'2-9'9'9'9'2-2'2'2'2-9'9-9'9'9'2'2- 

ooooooooooooooooooooooooooooooooooooooooooooo 


%Data  import  code  for  muon  fission  experiment 
%Written  by  Lt  Logan  Brandt 
%20  November  2014 


clear; 

clc; 


2'2'2-2'2'2'9'9'9-9'9'2-2'2-2'2-2'2'9'9'9'9'2-2' 

oooooooooooooooooooooooo 


%  Initialize  Variables  % 


9'2'9-2'2-2'9-9'2'2'2-2'2'2'2'2-2'2-9'2'2'9--2-2' 

oooooooooooooooooooooooo 


Nal_channels  =  4; 
neutron.channels  =  8; 

Fe_runs  =  6; 

PbLruns  =  3; 
Background_runs  =  2; 
Pb3_runs  =  5; 


9'2-9'2'9-2'9'9'9'9'9'9'2-2'9'9'2'9'9'9'9'9'2'2-2'2'9'9-9'2'9'9'9'9'9-9'2'9'9'9'9-9'9'9'2'2-9'2'9'9'9'9'9'9'2-9'2'9'9'9'9-9'9-9'9'2-2-2-2'2'9'9'9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%  IMPORTANT: 


Run  each  section 


individually 


and  save  the  workspace 


as 
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%  *type_of_material*_Data .mat .  This  file  type  is  a  required  input  for 
the 

%  analysis  code. 


2'9-2'9-2'2-9'9-9'9'9'2'9'2-2-2'9-9'2-9'9--9'2'9'2'9-2'2-9'2-9'9'9'2'9'2'9-2'2-2'2-9'9'9'2'9-2'9-2'2-2'2-9'9'9'2'9'2'9-2'2-9'9'9'9--9- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


9'2'9-2'9-2'2'9'2'2'9-9'2'2'2'9-2'2'9'9'2'2'2-2'9-2'9-9' 

oooooooooooooooooooooooooooo 


%  Import  Iron  Fission  Data  % 


9'2'2-2'2'2'9-9'9-9'9'2-2'2-2'2'2'2'2'9'9'9'2-2'2-2'2'9- 

oooooooooooooooooooooooooooo 


for  i  =  1:  NaLchannels 
Nal  =  [] ; 
for  j  =  3:Fe_runs 


hold_file  =  load ( [ 1 /Users/LJBrandt/Documents/ThesisData/Fe/ 
Fe_PHA_0  0  ’  .  .  . 

num2str(j)  1 _ls_ 1  num2str(i-l)  ' . dat ' ] ) ; 

Nal  =  cat(l,  Nal,  hold.file) ; 


end 

NaI_Data{i}  =  Nal; 

end 


for  i 


0:  neutron_channels-l 


PSDO  =  [ ] ; 

PSD1  =  [] ; 

PSD2  =  [ ] ; 

for  j  =  3:  Fe_runs 
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holcLfileO  =  load ( [ ' /Users/LJBrandt/Documents/ThesisData/Fe/ 
Fe_PSD0_00  '  .  .  . 

num2str(j)  ' _ls_'  num2str(i)  ' . dat ' ] ) ; 
hold_filel  =  load ( [ ' /Users/LJBrandt/Documents/ThesisData/Fe/ 
Fe_PSDl_00  '  .  .  . 

num2str(j)  ' _ls_'  num2str(i)  ' . dat " ] ) ; 
hold_file2  =  load ( [ ' /Users/LJBrandt/Documents/ThesisData/Fe/ 
Fe_P  SD2  _0  0  '  .  .  . 

num2str(j)  ' _ls_'  num2str(i)  ' . dat  ] ) ; 


PSDO 

=  cat  ( 1 , 

PSDO, 

hold_f ileO ) ; 

PSD1 

=  cat  ( 1 , 

PSD1 , 

hold_f ilel ) ; 

PSD2 

=  cat  ( 1 , 

PSD2 , 

hold_f ile2 ) ; 

end 

PSDO_Data{i+l}  =  PSDO; 

PSDl_Data{i+l}  =  PSD1; 

PSD2_Data{i+l}  =  PSD2 ; 

end 

clear  Background_runs  Fe_runs  hold_file  hold_fileO  hold_filel  hold_file2 

i  j  Nal  Nal_channels  neutronxhannels  Pbl_runs  Pb3_runs  PSDO  PSD1 
PSD2 

o,  o_ 
o  o 


9'2'9-2'9-2'9-9'9'9'9'9'2'9'2'9-2'9-9'9'9'9'9-2'9-2'9-9'9-9'2'9'9'9-2'9- 

oooooooooooooooooooooooooooooooooooo 


%  Import  1  Block  Lead  Fission  Data  % 


9'2'2-2'2'2'9-9'9'9'9'2-2'2-2'2'2'2'9'9'9'9'2-2'2-2'2'9'9'9'9'9'9'2-2'2- 

oooooooooooooooooooooooooooooooooooo 
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for  i 


1:  NaLchannels 


Nal  =  [] ; 

for  j  =  1:  Pbl_runs 


hold-file  =  load ( [ ' /Users/LJBrandt/Documents/ThesisData/Pbl/ 
Pbl_PHA_0  0  '  .  .  . 

num2str(j)  '_ls_'  num2str(i-l)  ' . dat ' ] ) ; 

Nal  =  cat(l,  Nal,  hold.file) ; 


end 

NaI_Data{i}  =  Nal; 

end 

for  i  =  0:  neutron_channels-l 
PSDO  =  [ ] ; 

PSD1  =  [] ; 

PSD2  =  [ ] ; 

for  j  =  1:  Pbbruns 

hold_file0  =  load ( [ ' /Users/LJBrandt/Documents/ThesisData/Pbl/ 
Pbl_PSD0_0  0  '  .  .  . 

num2str(j)  ' _ls_ 1  num2str(i)  '.dat"]); 
hold_filel  =  load ( [ ' /Users/LJBrandt/Documents/ThesisData/Pbl/ 
Pbl_PSDl_0  0  '  .  .  . 

num2str(j)  '_ls_'  num2str(i)  '.dat"]); 
hold_file2  =  load ( [ ' /Users/LJBrandt/Documents/ThesisData/Pbl/ 
Pbl_PSD2_0  0  '  .  .  . 

num2str(j)  ' _ls_ '  num2str(i)  '.dat  ]); 


PSDO 

=  cat  ( 1 , 

PSDO, 

hold_f ileO ) ; 

PSD1 

=  cat  ( 1 , 

PSD1 , 

hold-filel) ; 

PSD2 

=  cat  ( 1 , 

PSD2 , 

hold_f ile2 ) ; 
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end 


PSDO_Data{i+l}  =  PSDO; 

PSDl_Data{i+l}  =  PSD1; 

PSD2_Data{i+l}  =  PSD2 ; 

end 

clear  Background_runs  Fe_runs  hold-file  hold_fileO  hold_filel  hold_file2 

i  j  Nal  Nal_channels  neutronxhannels  Pbl_runs  Pb3_runs  PSDO  PSD1 
PSD2 

o,  o, 
o  o 


9'2'2-9'2'2-9-9'9'9'9'2-2'2-2'2'9'9'9'9'9'9'2-2'2-9- 

oooooooooooooooooooooooooo 


%  Import  Background  Data  % 


9'2'9'2'2-2'9-9'9-9'9'9'2'9-2'9-2'9-9'9'9'9'9-2'9-2' 

oooooooooooooooooooooooooo 


for  i  =  1:  NaLchannels 
Nal  =  [] ; 

for  j  =  1:  Background_runs 


hold-file  =  load  (  [ ' /Users/LJBrandt/Documents/ThesisData/ 
Background/Background_P HA_00  '  .  .  . 
num2str(j)  ' _ls_ 1  num2str(i-l)  ' . dat ' ] ) ; 

Nal  =  cat(l,  Nal,  hold-file) ; 


end 

NaI_Data{i}  =  Nal; 

end 

for  i  =  0:  neutron_channels-l 
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PSDO  =  [ ] ; 


PSD1  =  [] ; 

PSD2  =  [ ] ; 

for  j  =  1:  Background-runs 

holcLfileO  =  load  (  [ ' /Users/LJBrandt/Documents/ThesisData/ 
Background/Background_PSDO_0 0  '  .  .  . 
num2str(j)  ' _ls_ 1  num2str(i)  ' . dat  ]); 
hold_filel  =  load ( [ ' /Users/LJBrandt/Documents/ThesisData/ 
Background/Background_PSDl_00 ' . . . 
num2str(j)  '_ls_'  num2str(i)  '.dat  ]); 
hold_file2  =  load ( [ ' /Users/LJBrandt/Documents/ThesisData/ 
Background/Background_PSD2_0 0 ' . . . 


num2str ( j ; 

i  '_ls. 

_ '  num2str (i) 

PSDO 

=  cat  ( 1 , 

PSDO, 

hold_f ileO ) ; 

PSD1 

=  cat  ( 1 , 

PSD1 , 

hold-filel) ; 

PSD2 

=  cat  ( 1 , 

PSD2 , 

hold_f ile2 ) ; 

end 

PSDO_Data{i+l}  =  PSDO; 

PSDl_Data{i+l}  =  PSD1; 

PSD2_Data{i+l}  =  PSD2 ; 

end 

clear  Background_runs  Fe_runs  hold-file  hold_fileO  hold_filel  hold_file2 


i  j  Nal  Nal_channels  neutronxhannels  Pbl_runs  Pb3_runs  PSDO  PSD1 
PSD2 

o,  o. 
o  o 


9'2'9-9'2'2'9'9'9-9'9'2-2'2-2'2'9'9'9'9'9'9'2-2-2-9'2'9'9-9'9'9'9'2-2'2- 

oooooooooooooooooooooooooooooooooooo 


%  Import  3  Block  Lead  Fission  Data  % 


56 


for  i 


1:  Nal_channels 


Nal  =  [] ; 

for  j  =  1:  Pb3_runs 

hold-file  =  load ( [ ' /Users/LJBrandt/Documents/ThesisData/Pb3/ 
Pb3_PHA_0  0  '  .  .  . 

num2str(j)  ' _ls_ 1  num2str(i-l)  ' . dat ' ] ) ; 

Nal  =  cat(l,  Nal,  hold.file) ; 

end 

NaI_Data{ i}  =  Nal; 

end 

for  i  =  0:  neutron_channels-l 
PSDO  =  [ ] ; 

PSD1  =  [] ; 

PSD2  =  [ ] ; 

for  j  =  1:  Pb3_runs 

hold_file0  =  load ( [ ' /Users/LJBrandt/Documents/ThesisData/Pb3/ 
Pb3_PSD0_0  0  '  .  .  . 

num2str(j)  '_ls_'  num2str(i)  '.dat  ]); 
hold_filel  =  load ( [ ' /Users/LJBrandt/Documents/ThesisData/Pb3/ 
Pb3_PSDl_0  0  '  .  .  . 

num2str(j)  '_ls_'  num2str(i)  '.dat  ]); 
hold_file2  =  load ( [ ' /Users/LJBrandt/Documents/ThesisData/Pb3/ 
Pb3_PSD2_0  0  '  .  .  . 

num2str(j)  '  dsd  num2str(i)  '.dat  ]); 
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PSDO 

=  cat  ( 1 , 

PSDO, 

hold_f ileO ) ; 

PSD1 

=  cat  ( 1 , 

PSD1 , 

hold_filel) ; 

PSD2 

=  cat  ( 1 , 

PSD2 , 

hold_f ile2 ) ; 

end 

PSDO_Data{i+l}  =  PSDO; 

PSDl_Data{i+l}  =  PSD1; 

PSD2_Data{i+l}  =  PSD2 ; 

end 

clear  Background-runs  Fe.runs  hold-file  hold_fileO  hold_filel  hold_file2 

i  j  Nal  Nal_channels  neutron.channels  Pbl.runs  Pb3_runs  PSDO  PSD1 
PSD2 


9'2'2-2'2'9'9-9'9-9'9'9-2'2-2'9'9'9'9'9'9'9'2-2'2-9'2'9'9-9'9'9'9'2-2'2-2'9'9'9'9'9'9'2'2- 

ooooooooooooooooooooooooooooooooooooooooooooo 


9'2'2-2'2'2'9'2'9-2'9'2-2'2-2'2'2-2'9'9'9'9'2-2'2-2'2'9'2'2'9'9'9'2-2'2'2'2-2'9-9'9'9'2'2- 

ooooooooooooooooooooooooooooooooooooooooooooo 


%Bethe  Block  Calculations 


2'2'2-2'9-2'2-9'2'9'9'9-2'2-2'2'2'2'9'2-2'9'2-2'9-2'9-9'9-9'9'9'2'9-2'9-2'2-9'2-9'9--9'2'9- 

ooooooooooooooooooooooooooooooooooooooooooooo 


9'2'9'9'9-2'9-9'9-9'9'9-2'9-2'9-2'9-9'9'9'9'9-2'9'2'9-9'9-9'9'9'9'9-2'2'2'9-9'9-9'9'9'2'9- 

ooooooooooooooooooooooooooooooooooooooooooooo 


%Written  by  Maj  Van  Dyk 
%Edited  by  Lt  Brandt 


%General  Parameters 
c  =  2 . 997 92458*10 " 8;  %meters/sec 

re  =  2 . 817 94 0325*1 0 ~  (-13 )  ;  %radius  of  an  electron  in  centimeters 
me  =  0 . 510998918/c“2; %mass  of  an  electron  in  MeV 
Mpart  =  105 . 658372 /c~2 ;  %mass  of  a  muon  in  MeV 
Na  =  6 . 0221415*10~23;  %Avagadros  number 
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zpart  =  1 .  ; 


%NaI  Parameters 
ZdivA  =  0.42697; 

rho  =  3 . 6 67 ;  %density  of  Nal=3.667  g/cm~3 

Ionization  =  452  .  *10"'  (-6)  ;  %MeV 

thickness  =  5.08;%cm 

cbar  =  6.0572; 

a  =  0.12516; 

xl  =  3.5920; 

smallk  =  3.0398; 

%Stainless  Steel  Parameters 
ZdivA_SS  =  0.46556; 
rho_SS  =  7.8740; 

Ionization_SS  =  28 6 . *10 “  (-6) ;  %MeV 
thickness_SS_2  =  0.1016*2;%cm 
thickness_SS_l  =  0.1016;%cm 
cbar.SS  =  4.2911; 
a_SS  =  0.14680; 
x 1 _S  S  =  3.1531; 
smallk_SS  =  2.9632; 

%Teflon  Parameters 
ZdivA_Teflon  =  0.47992; 
rho_Teflon  =  2.2; 

Ionization_Tef Ion  =99 . 1*10“  (-6) ;  %MeV 

thickness_Tef Ion  =  0.2159;%cm 

cbar_Teflon  =  3.4161; 

a_Teflon  =  0.10606; 

xl_Teflon  =  2.7404; 

smallk_Tef Ion  =  3.4046; 
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dElBB.Nal 


o, 
o 

%  dE2BB_NaI 
%  dE3BB_NaI 
%  dE4BB_NaI 

incrl  =  1; 

for  KEnergyl_iterate=140 : 1 : 5000 

%Find  Emu2BB 


%Find  decrease  in  energy  due  to  Nal  1 
muon_energy=KEnergyl  .iterate ; 

vel=sqrt (-c~2* ( (1/ (muon.energy/ (Mpart*c"2) +1) ~2) -1) ) ; 

gamma  =  l/(sqrt(l  -  (vel/c) " 2 ) ) ; 
beta  =  vel/c; 

Tmax  =  (2*me*c''2*beta''2*gamma''2)  /  (1  +  (2*gamma*me/Mpart) 
+  (me/Mpart)  “2) ; 
momen  =  Mpart*beta*gamma*c; 
frac  =  momen/ (Mpart*c) ; 
x  =  loglO (frac) ; 

densitycorrection  =  2*log(10)*x  -  cbar  +  a*  (xl  -  x) 
smallk ; 

K  =  4*pi*Na*re~2*me*c,'2; 

dElBB.Nal (incrl )  =  K*zpart~2* (ZdivA) * (l/beta~2) * ( (1/2) * 
log  (  (2*me*c''2*beta~2*gamma''2*Tmax)  /Ionization-^  )  - 
beta-'2  -  densitycorrection/2 )  *thickness*rho;  %MeV 

KEnergy l_BB_NaI=KEnergyl_iterate-dElBB_NaI (incrl) ; 
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%Find  decrease  in  energy  due  to  Teflon 

muon_energy=KEnergy  l_BB_NaI ; 

vel=sqrt  (-c"2*  (  (1/  (muon.energy/  (Mpart*c'2)  +1)  ~2)  -1)  )  ; 

gamma  =  l/(sqrt(l  -  (vel/c) " 2 ) ) ; 

beta  =  vel/c; 

Tmax  =  (2*me*c~2*beta''2*gamma''2)  /  (1  +  (2*gamma*me/Mpart) 
+  (me/Mpart) ~2) ; 

momen  =  Mpart*beta*gamma*c; 

frac  =  momen/ (Mpart*c) ; 

x  =  loglO (frac) ; 

densitycorrection  =  2*log(10)*x  -  cbar.Teflon  +  a.Teflon 
* (xl.Tef Ion  -  x) smallk.Tef Ion; 

K  =  4*pi*Na*re~2*me*c~2; 

dElBB_Tef londower  =  K*zpart~2*  (ZdivA_Tef Ion)  *  (l/beta''2) 
*  (  (1/2)  *log  (  (2*me*c''2*beta''2*gamma''2*Tmax)  / 
Ionization_Teflon"'2)  -  beta''2  -  densitycorrection/2 ) 
*thi ckne ss _Tefl on *rho .Teflon;  %MeV 

KEnergyl_BB_Teflon_lower=KEnergyl_BB_NaI- 
dElBB_Tef  Ion  .lower ; 


%Find  decrease  in  energy  due  to  Stainless  Steel 
muon_energy=KEnergy  l_BB_Tef  londower  ; 

vel=sqrt (-c"2* ( (1/ (muon.energy/ (Mpart*c"2) +1) ~2) -1) )  ; 

gamma  =  l/(sqrt(l  -  (vel/c) ~ 2 ) ) ; 
beta  =  vel/c; 

Tmax  =  (2*me*c''2*beta''2*gamma''2)  /  (1  +  (2*gamma*me/Mpart) 
+  (me/Mpart) ~2) ; 
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momen  =  Mpart*beta*gamma*c; 
frac  =  momen/ (Mpart*c) ; 
x  =  loglO (frac) ; 

densitycorrection  =  2*log(10)*x  -  cbar_SS  +  a_SS*(xl_SS 
-  x)  smallk_SS ; 

K  =  4*pi*Na*re~2*me*c~2; 

dElBB_SS  =  K*zpart~2* (ZdivA_SS) * (l/beta'2) * ( (1/2) *log 
(  (2*me*c''2*beta~2*gamma~2*Tmax)  /Ionization_SS~2)  - 
beta~2  -  densitycorrection/2 ) *thickness_SS_2*rho_SS; 

%MeV 

KEnergy  l_BB_SS=KEnergy  l_BB_Tef  lon_lower-dElBB_SS  ; 

%Find  decrease  in  energy  due  to  Teflon 
muon_energy=KEnergyl_BB_SS  ; 

vel=sqrt (-c"2* ( (1/ (muon_energy/ (Mpart*c"2) +1) ~2) -1) ) ; 

gamma  =  l/(sqrt(l  -  (vel/c) " 2 ) )  ; 
beta  =  vel/c; 

Tmax  =  (2*me*c''2*beta''2*gamma''2)  /  (1  +  (2*gamma*me/Mpart) 
+  (me/Mpart)  “2) ; 
momen  =  Mpart*beta*gamma*c; 
frac  =  momen/ (Mpart*c) ; 
x  =  loglO (frac) ; 

densitycorrection  =  2*log(10)*x  -  cbar_Teflon  +  a_Teflon 
* (xl_Tef Ion  -  x) ~smallk_Tef Ion; 

K  =  4*pi*Na*re~2*me*c''2; 

dE2BB_Tef lon.upper  =  K*zpart~2*  (ZdivA_Tef Ion)  *  (l/beta''2) 
*  (  (1/2)  *log  (  (2*me*c''2*beta''2*gamma''2*Tmax)  / 
Ionization_Teflon''2)  -  beta'2  -  densitycorrection/2) 
*thi ckne ss _Tefl on *rho .Teflon;  %MeV 
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KEnergy2_BB=KEnergyl_BB_SS-dE2BB_Tef  lon.upper  ; 


find  Emu3BB 


%Find  decrease  in  energy  due  to  Nal  2 
muon_energy=KEnergy2_BB; 

vel=sqrt  (-c~2*  (  (1/  (muon.energy/  (Mpart*c'2)  +1)  ~2)  -1)  )  ; 

gamma  =  l/(sqrt(l  -  (vel/c) ~ 2 ) ) ; 
beta  =  vel/c; 

Tmax  =  (2*me*c~2*beta~2*gamma~2)  /  (1  +  (2*gamma*me/Mpart) 
+  (me/Mpart) ~2) ; 
momen  =  Mpart*beta*gamma*c; 
frac  =  momen/ (Mpart*c) ; 
x  =  loglO (frac) ; 

densitycorrection  =  2*log(10)*x  -  cbar  +  a* (xl  -  x) 
smallk ; 

K  =  4*pi*Na*re''2*me*c"2; 

dE2BB_NaI (incrl )  =  K*zpart~2* (ZdivA) * (l/beta~2) * ( (1/2) * 
log  (  (2*me*c,'2*beta''2*gamma''2*Tmax)  /Ionization~2  )  - 
beta"'2  -  densitycorrection/2 )  *thickness*rho;  %MeV 

KEnergy2_BB_NaI=KEnergy2_BB-dE2BB_NaI  (incrl)  ; 


%Find  decrease  in  energy  due  to  Teflon 
muon_energy=KEnergy2_BB_NaI ; 

vel=sqrt  (-c"2*  (  (1/  (muon_energy/  (Mpart*c,'2)  +1)  ~2)  -1)  )  ; 
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gamma  =  l/(sqrt(l  -  (vel/c) ~ 2 ) ) ; 
beta  =  vel/c; 

Tmax  =  (2*me*c~2*beta~2*gamma~2)  /  (1  +  (2*gamma*me/Mpart) 
+  (me/Mpart) ~2) ; 
momen  =  Mpart*beta*gamma*c; 
frac  =  momen/ (Mpart*c) ; 
x  =  loglO (frac) ; 

densitycorrection  =  2*log(10)*x  -  cbar_Teflon  +  a_Teflon 
* (xl_Tef Ion  -  x) ~ smallk_Tef Ion; 

K  =  4*pi*Na*re''2*me*c"2; 

dE2BB_Tef lon_lower  =  K*zpart~2*  (ZdivA_Tef Ion)  *  (l/beta"2) 
*  (  (1/2)  *log  (  (2*me*c''2*beta''2*gamma''2*Tmax)  / 
Ionization_Teflon"'2)  -  beta''2  -  densitycorrection/2 ) 
*thi ckne ss _Tefl on *rho .Teflon;  %MeV 

KEnergy2_BB_Teflon_lower=KEnergy2_BB_NaI- 
dE2BB_Tef  londower ; 


%Find  decrease  in  energy  due  to  Stainless  Steel 
muon_energy=KEnergy2_BB_Tef  lon.lower  ; 

vel=sqrt  (-c~2*  (  (1/  (muon_energy/  (Mpart*c,'2)  +1)  ~2)  -1)  )  ; 

gamma  =  l/(sqrt(l  -  (vel/c) ~ 2 ) ) ; 
beta  =  vel/c; 

Tmax  =  (2*me*c''2*beta,'2*gamma''2)  /  (1  +  (2*gamma*me/Mpart) 
+  (me/Mpart) ~2) ; 
momen  =  Mpart*beta*gamma*c; 
frac  =  momen/ (Mpart*c) ; 
x  =  loglO (frac) ; 

densitycorrection  =  2*log(10)*x  -  cbar_SS  +  a_SS*(xl_SS 
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-  x)  smal lk_SS ; 


K  =  4*pi*Na*re''2*me*c"2; 

dE2BB_SS  =  K*zpart~2* (ZdivA_SS) * (l/beta'2) * ( (1/2) *log 
(  (2*me*c''2*beta~2*gamma~2*Tmax)  /Ionization_SS~2)  - 
beta~2  -  densitycorrection/2 ) *thickness_SS_2*rho_SS; 
%MeV 

KEnergy2_BB_SS=KEnergy2_BB_Tef  lon_lower-dE2BB_SS  ; 


%Find  decrease  in  energy  due  to  Teflon 

muon_energy=KEnergy2_BB_SS  ; 

vel=sqrt (-c~2* ( (1/ (muon.energy/ (Mpart*c"2) +1) ~2) -1) ) ; 

gamma  =  l/(sqrt(l  -  (vel/c) "2) ) ; 

beta  =  vel/c; 

Tmax  =  (2*me*c~2*beta''2*gamma~2)  /  (1  +  (2*gamma*me/Mpart) 
+  (me/Mpart) ~2) ; 

momen  =  Mpart*beta*gamma*c; 

frac  =  momen/ (Mpart*c) ; 

x  =  loglO (frac) ; 

densitycorrection  =  2*log(10)*x  -  cbar.Teflon  +  a.Teflon 
* (xl_Tef Ion  -  x) " smallk_Tef Ion; 

K  =  4*pi*Na*re~2*me*c~2; 

dE3BB_Tef lonmpper  =  K*zpart~2*  (ZdivA_Tef Ion)  *  (l/beta"2) 
*  (  (1/2)  *log  (  (2*me*c~2*beta~2*gamma~2*Tmax)  / 
Ionization_Teflon"2)  -  beta''2  -  densitycorrection/2) 
*thi ckne ss _Tefl on *rho .Teflon;  %MeV 

KEnergy3_BB=KEnergy2_BB_SS-dE3BB_Tef  lonmpper  ; 
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find  Emu4BB 


%Find  decrease  in  energy  due  to  Nal  3 
muon_energy=  KEnergy3_BB; 

vel=sqrt  (-c~2*  (  (1/  (muon_energy/  (Mpart*c'2)  +1)  ~2)  -1)  )  ; 

gamma  =  l/(sqrt(l  -  (vel/c) ~ 2 ) ) ; 
beta  =  vel/c; 

Tmax  =  (2*me*c~2*beta~2*gamma~2) / (1  +  (2*gamma*me/Mpart) 
+  (me/Mpart) ~2) ; 
momen  =  Mpart*beta*gamma*c; 
frac  =  momen/ (Mpart*c) ; 
x  =  loglO (frac) ; 

densitycorrection  =  2*log(10)*x  -  cbar  +  a*  (xl  -  x) ~ 
smallk; 

K  =  4*pi*Na*re''2*me*c"2; 

dE3BB_NaI  (incrl )  =  K*zpart"'2*  (ZdivA)  *  (l/beta"'2)  *  (  (1/2)  * 
log  (  (2*me*c,'2*beta''2*gamma~2*Tmax)  /Ionization~2  )  - 
beta"'2  -  densitycorrection/2 )  *thickness*rho;  %MeV 

KEnergy 3_BB_NaI=KEnergy3_BB-dE3BB_NaI  (incrl)  ; 


%Find  decrease  in  energy  due  to  Teflon 
muon_energy=KEnergy3_BB_NaI ; 

vel=sqrt (-c"2* ( (1/ (muon.energy/ (Mpart*c"2) +1) ~2) -1) ) ; 

gamma  =  l/(sqrt(l  -  (vel/c) " 2 ) ) ; 
beta  =  vel/c; 

Tmax  =  (2*me*c''2*beta''2*gamma''2)  /  (1  +  (2*gamma*me/Mpart) 
+  (me/Mpart)  “2) ; 
momen  =  Mpart*beta*gamma*c; 
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frac  =  momen/ (Mpart*c) ; 

x  =  loglO (frac) ; 

densitycorrection  =  2*log(10)*x  -  cbar.Teflon  +  a.Teflon 
* (xl_Tef Ion  -  x) ~ smallk_Tef Ion; 

K  =  4*pi*Na*re~2*me*c~2; 

dE3BB_Tef lon.lower  =  K*zpart^2*  (ZdivA_Tef Ion)  *  (l/beta"2) 
*  (  (1/2)  *log  (  (2*me*c~2*beta~2*gamma~2*Tmax)  / 
Ionization_Teflon"2)  -  beta'2  -  densitycorrection/2 ) 
*thi ckne ss _Tefl on *rho -Teflon;  %MeV 

KEnergy3_BB_Teflon_lower=KEnergy3_BB_NaI- 
dE3BB_Tef Ion .lower ; 


%  Find  decrease  in  energy  due  to  Stainless  Steel 
muon_energy=KEnergy3_BB_Tef  londower  ; 

vel=sqrt  (-c~2*  (  (1/  (muon_energy/  (Mpart*c"2)  +1)  ~2)  -1)  )  ; 

gamma  =  l/(sqrt(l  -  (vel/c) "2) ) ; 
beta  =  vel/c; 

Tmax  =  (2*me*c~2*beta''2*gamma''2)  /  (1  +  (2*gamma*me/Mpart) 
+  (me/Mpart)  “2) ; 
momen  =  Mpart*beta*gamma*c; 
frac  =  momen/ (Mpart*c) ; 
x  =  loglO (frac) ; 

densitycorrection  =  2*log(10)*x  -  cbar_SS  +  a_SS*(xl_SS 
-  x)  smallk_SS ; 

K  =  4*pi*Na*re~2*me*c~2; 

dE3BB_SS  =  K*zpart~2* (ZdivA_SS) * (l/beta'2) * ( (1/2) *log 
(  (2*me*c''2*beta~2*gamma~2*Tmax)  /Ionization_SS~2)  - 
beta~2  -  densitycorrection/2  )  *thicknes  s_SS_2  *  rho_SS  ; 

%MeV 
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KEnergy3_BB_SS=KEnergy3_BB_Tef  lon_lower-dE3BB_SS  ; 


%Find  decrease  in  energy  due  to  Teflon 

muon_energy=KEnergy3_BB_SS  ; 

vel=sqrt  (-c"2*  (  (1/  (muon.energy/  (Mpart*c"2)  +1)  ~2)  -1)  )  ; 

gamma  =  l/(sqrt(l  -  (vel/c) ~ 2 ) ) ; 

beta  =  vel/c; 

Tmax  =  (2*me*c~2*beta~2*gamma''2)  /  (1  +  (2*gamma*me/Mpart) 
+  (me/Mpart)  “2) ; 

momen  =  Mpart*beta*gamma*c; 

frac  =  momen/ (Mpart*c) ; 

x  =  loglO (frac) ; 

densitycorrection  =  2*log(10)*x  -  cbar.Teflon  +  a.Teflon 
* (xl_Tef Ion  -  x) ' smallk_Tef Ion; 

K  =  4*pi*Na*re''2*me*c''2; 

dE4BB_Tef lon_upper  =  K*zpart,'2*  (ZdivA_Tef Ion)  *  (l/beta''2) 
*  (  (1/2)  *log  (  (2*me*c''2*beta''2*gamma''2*Tmax)  / 
Ionization_Teflon''2)  -  beta'2  -  densitycorrection/2 ) 
*thi ckne ss _Tefl on *rho .Teflon;  %MeV 

KEnergy4_BB=KEnergy3_BB_SS-dE4BB_Tef  lonmpper  ; 


find  muon  exiting  energy 

%Find  decrease  in  energy  due  to  Nal  4 
muon_energy=  KEnergy4_BB; 

vel=sqrt (-c"2* ( (1/ (muon.energy/ (Mpart*c'2) +1) ~2) -1) ) ; 

gamma  =  l/(sqrt(l  -  (vel/c) ~ 2 ) ) ; 
beta  =  vel/c; 
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Tmax  =  (2*me*c~2*beta~2*gamma~2)  /  (1  +  (2*gamma*me/Mpart) 
+  (me/Mpart) ~2) ; 

momen  =  Mpart*beta*gamma*c; 

frac  =  momen/ (Mpart*c) ; 

x  =  loglO (frac) ; 

densitycorrection  =  2*log(10)*x  -  cbar  +  a*  (xl  -  x) ~ 
smallk ; 

K  =  4*pi*Na*re~2*me*c"2; 

dE4BB_NaI (incrl )  =  K*zpart~2* (ZdivA) * (l/beta~2) * ( (1/2) * 
log  (  (2*me*c,'2*beta''2*gamma~2*Tmax)  /Ionization~2  )  - 
beta~2  -  densitycorrection/2 ) *thickness*rho;  %MeV 

KEnergy 4_BB_NaI=KEnergy4_BB-dE4BB_NaI  (incrl)  ; 


%Find  decrease  in  energy  due  to  Teflon 
muon_energy=KEnergy4_BB_NaI ; 

vel=sqrt  (-c"2*  (  (1/  (muon_energy/  (Mpart*c,'2)  +1)  ~2)  -1)  )  ; 

gamma  =  l/(sqrt(l  -  (vel/c) ”2 ) ) ; 
beta  =  vel/c; 

Tmax  =  (2*me*c''2*beta''2*gamma''2)  /  (1  +  (2*gamma*me/Mpart) 
+  (me/Mpart) ~2) ; 
momen  =  Mpart*beta*gamma*c; 
frac  =  momen/ (Mpart*c) ; 
x  =  loglO (frac) ; 

densitycorrection  =  2*log(10)*x  -  cbar_Teflon  +  a_Teflon 
* (xl_Tef Ion  -  x) " smallk_Tef Ion; 

K  =  4*pi*Na*re''2*me*c"2; 

dE4BB_Tef londower  =  K*zpart'2*  (ZdivA_Tef Ion)  *  (l/beta''2) 
*  (  (1/2)  *log  (  (2*me*c''2*beta''2*gamma''2*Tmax)  / 
Ionization_Teflon"'2)  -  beta''2  -  densitycorrection/2 ) 
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*thi ckne ss _Tefl on *rho .Teflon;  %MeV 


KEnergy4_BB_Tef  lon_lower=KEnergy4_BB_NaI  - 
dE4BB_Tef Ion .lower ; 


%  Find  decrease  in  energy  due  to  Stainless  Steel 
muon_energy=KEnergy4_BB_Tef  lon.lower  ; 

vel=sqrt (-c"2* ( (1/ (muon_energy/ (Mpart*c"2) +1) ~2) -1) )  ; 

gamma  =  l/(sqrt(l  -  (vel/c) ~ 2 ) ) ; 
beta  =  vel/c; 

Tmax  =  (2*me*c~2*beta''2*gamma''2)  /  (1  +  (2*gamma*me/Mpart) 
+  (me/Mpart) ~2) ; 
momen  =  Mpart*beta*gamma*c; 
frac  =  momen/ (Mpart*c) ; 
x  =  loglO (frac) ; 

densitycorrection  =  2*log(10)*x  -  cbar_SS  +  a_SS*(xl_SS 
-  x)  smallk_SS ; 

K  =  4*pi*Na*re''2*me*c"2; 

dE4BB_SS  =  K*zpart"2* (ZdivA_SS) * (l/beta'2) * ( (1/2) *log 
(  (2*me*c''2*beta~2*gamma~2*Tmax)  /Ionization_SS-'2)  - 
beta"'2  -  densitycorrection/2)  *thickness_SS_l*rho_SS; 
%MeV 

KEnergy_exit_BB  (incrl)  =KEnergy4_BB_Tef lon_lower“dE4BB_SS 


incrl=incrl+l; 


end 
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Save  workspace  as  BetheBlock .mat 


9'2'9'2'2-2'9-9'9-9'9'9'2'9'2'9-2'9-9'9'9'9'9-2'9-2'9-9'9-9'9'9'9' 

ooooooooooooooooooooooooooooooooo 


2'2'2-2'2'2'9'9'9-2'9'2-2'2-2'2'2'2'9'9'9'9'2-2'2-2'2'9'2'9'9'2'9'2-2'2'2- 

ooooooooooooooooooooooooooooooooooooo 


9'2'2-9'2'9'9-9'9-9'9'9-2'2-2'2'2'9'9'9'9'9'2-2'2-9'2'9-9'9'9'9'9'2-9'2'2' 

ooooooooooooooooooooooooooooooooooooo 


%MAIN  PORTION  OF  CODE 


%Performs  Neutron  Timing  calculations 
%and  Nal  muon  energy  determination 


9'2'2-2'2'2'9'9'9-2'9'2-2'2-2'2'2'2'9'9'9'9'2-2'2-2'2'9'2'2'9'9'9'2-2'2'2- 

ooooooooooooooooooooooooooooooooooooo 


9'2'2-2'2'2'9'9'9-9'9'2-2'2-2'2'2'2'9'9'9'9'2-2'2-2'2'9'9-2-9'9'9'2-2'2'2- 

ooooooooooooooooooooooooooooooooooooo 


%Written  by  Maj  Greg  Van  Dyk 
%07  November  2013 
%Edited  by  Lt  Logan  Brandt 
%November  2014 


close 

all 

clear 

all 

clc 

9'2'2-2'2'2'9'9'9-9'9'9-9'2-2'9'9'9'9'9'9'9'2-2'2-9-2'9'9-9'9'9'9'2-9'2'2' 

ooooooooooooooooooooooooooooooooooooo 


%  File  Characteristics  for  input  data 


9'2'9-9'9-2'9-9'9-9'9'9-2'9-2'9-2'9-9'9'9'9'9-2'9-2'9-9'9-9'2'9'9'9-2'2'2' 

ooooooooooooooooooooooooooooooooooooo 


total _files_NaI=4; 
total_f iles_neutrons=8; 
concident_acquisition_window=5000 ; 
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%Coincident  window  for  neutron 
%in  nanoseconds  dividied  by  4; 


digitizers,  take 
20  microsec=5000 


the  time 

50  microsec=12500 


%Initializations 
counter=l ; 
counterl  =  l ; 
counter3=l ; 
counter5=l ; 
incr=l ; 
incrl=l ; 

energy4_exiting_hist= [ ] ; 

progress=0 ; 


9'2'2-2'2'9'9-9'9'9'9'2-2'2-2'9'9'9'9'9'9'9'2-2'2-9'2'9'9'9'9'9'9'2-2'2'2'2'9'9-9-9'9'9'2-2'2'2'2'9'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%Import  Data 


9'2'9'2'2'2'9-9'9-9'9'9-2'9-2'9-2'9-9'9'9'9'9-2'9-9'2'9'9-9'9'9'9'9-2'2'2'9-9'9-9'9'9'2'9-2'9-2'2'9'9-9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooo 


load ( 1 Background_Data . mat 1 ) 
material 


%load ( 

'  Fe_Data  .  mat  ' 

%load ( 

'  Pbl_Data .  mat 

%load ( 

'  Pb2_Data .  mat 

%load ( 

'  Pb3_Data .  mat 

%%%%  Change  with  new 


9'2'2-2'2'2'9-2'2'9'2'2-2'2-2'2'2'2'9'9'9'9'2-2'2-2'2'9'2'9'9'9'9'2-2'2'2-2'9'9-9'9'9'2-2-2'2'2'2'9'9-9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%  Data  Input  for  the  Neutron  Channels 


9'2'9'2'2'2'9-9'9-9'9'9'2'9'2'9-2'9-9'9'9'9'9'9'9'9'2'9'9-9'9'9'9'9-2'2'2'9-9'9-9'9'9'2'9-2'2-2'9-2'9-9'2' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%Import  the  first  neutron  digitizer  from  text  files  and  load  into 


matrices 


k=l  ; 
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for  i=2 : total_f iles.neutrons 


Data_neutronsl{i}=PSDO_Data{i}; 

if  numel (Data_neutronsl{i} )  >0 

PSD_neutrons 1 {k}=Data_neutrons 1 { i } ( : , 4) ; 
timing_neutronsl{k}=Data_neutronsl{i} ( : , 1) ; 
k=k+l; 

end 

end 

%Import  the  second  neutron  digitizer  from  text  files  and  load  into 
matrices 

k=l; 

for  i=2 : total_f iles.neutrons 

Data_neutrons2{i}=PSDl_Data{i}; 

if  numel (Data_neutrons2 { i } )  >0 

PSD_neutrons2{k}=Data_neutrons2{i} ( : , 4) ; 
timing_neutrons2{k}=Data_neutrons2{i} ( : , 1) ; 
k=k+l ; 

end 

end 

%Import  the  third  neutron  digitizer  from  text  files  and  load  into 
matrices 

k=l  ; 

for  i=2 : total_f iles.neutrons 

Data_neutrons3{i}=PSD2_Data{i}; 
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if  numel (Data_neutrons3{i})  >0 

PSD_neutrons3{k}=Data_neutrons3{i} ( : , 4) ; 
timing_neutrons3{k}=Data_neutrons3{i} ( : , 1) ; 
k=k+l ; 

end 

end 


9'2'9-2'9-2'9'9'9'9'9'9'2'9'9'9-2'9-9'9'9'9'9'2'9-2'9-9'9-9'9'9'2'9-2'9-2'9-9'9-9'9'9'2'9-2'9-2'9-9'9-9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%  Data  input  for  Nal(Tl)  detectors 


9'2'2-9'2'2'9'9'9-9'9'2-2'2-2'9'9'9'9'9'9'9'2-2'2-2'2'9'9'9'9'9'9'2-9'2'9'9'9'9'9'9'9'9'2-2'9-9'2'9'9-9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%Import  the  Nal  data  from  text  files  and  load  into  matrices 
for  i=0 : total_f iles_NaI-l 
j=i+l; 

Data{j}=NaI_Data{j}; 
pulse_height  s_NaI { j}=Data{j}(:,2) ; 
timing_Nal{ j}=Data{j} ( : , 1) ; 

end 


2'2'9-2'9-2'2-9'2'9'9-9'2'2'2'9-2'2-9'9'9'9'2-2'9-2'9-9'2'2'2'9'2'9-2'9-2'2'9'2-9'2'2'2'9-2'2-2'9-2'2-2'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%Data  Input  for  Frequency  Generator 


9'2'2-2'2'2'9-9'9-9'9'2'2'2-2'2'2'2'9'9'9'9'2-2'2-2'2'9'9'2'9'9'9'2-2'2'2'2'9'9-9'9'9'2-2-2'2'2'2'9'2-9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooo 


time_stampl=PSDO_Data{l} ( :  ,  1) ; 
check_sizel=numel (time_stampl ) ; 


time_stamp2=PSDl_Data{l} ( : , 1) ; 
check_size2=numel (time_stamp2 ) ; 
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time_stamp3=PSD2_Data{l} ( :  ,  1)  ; 
check_size3=numel (time_stamp3 ) ; 


display (' Data  Input  Complete') 


9'2'2-9'2'2'9'9'9-9'9'9-9'2-2'2'9'9'9'9'9'9'2-2'2-9'2'9'9-9'9'9'9'2-2'2'2'2'9'9-9'9'9'9'2-2'2'9'2'9'9'9'9'9'2'2-9'2'2'9'9'9-9'9'9'9'2-2'2- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%  Filtering  the  data  to  pull  out  just  one  single  event  per  time  stamp 
for 

%  frequency  generator  channels 

9'2'2-9'2'2'9'9'9'9'9'9'9'2-2'2'9'9'9'9'9'9'2-2'9-9'2'9'9'9'9'9'9'2-2'2'2'9'9'9-9'9'9'2'2-9'9'2'2'9'9'9'9'9'9'2-2'9'2'9'9'9-9'9'9'9'2-2'9- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


last_time_stampl=time_stampl (1) ; 
f iltered_time_stampsl  ( 1 )  =time_stampl  (1)  ; 
incr=2 ; 
counter=l ; 

for  m=2 : check_sizel 

current_time_stampl=time_stampl  (m)  ; 

timing_dif  fl  =  cur  rent  _time_st  amp  1 -last  _time_st  amp  1 ; 

if  abs (timing_dif f 1 )>concident_acquisition_window 

f  iltered_time_stampsl  (incr)  =current_t  ime_stampl ; 
incr=incr+l ; 

end 

last_time_stampl=time_stampl  (m)  ; 


end 
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Iast_time_stamp2=time_stamp2 (1) ; 
f iltered_time_stamps2  ( 1 )  =time_stamp2  (1)  ; 
incr=2 ; 

for  m=2 : check_size2 

current_time_stamp2=time_stamp2  (m)  ; 

timing_dif  f  2  =  cur rent _time_stamp2 -last _time_stamp2  ; 

if  abs (timing_dif  f2  )>concident_acquisition_window 

f  iltered_time_stamps2  (incr)  =current_t  ime_stamp2  ; 
incr=incr+l ; 

end 

Iast_time_stamp2=time_stamp2  (m)  ; 


end 


Iast_time_stamp3=time_stamp3  (1)  ; 
f iltered_time_stamps3  (1)  =time_stamp3  (1)  ; 
incr=2 ; 

for  m=2 : check_size3 

current_time_stamp3=time_stamp3  (m)  ; 

timing_dif  f  3=current_time_stamp3-last_time_stamp3  ; 

if  abs (timing_dif f3)>concident_acquisition_window 

f  i  Iter  ed_time_s  tamps  3  (incr)  =current_t  ime_stamp3  ; 
incr=incr+l ; 

end 

Iast_time_stamp3=time_stamp3  (m)  ; 
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end 


f iltered_t ime_stamps_t rans l=transpose  ( f ilteredrt ime_s tamps  1 )  ; 
f  iltered_t  ime_stamps_t  rans2=transpose  (  f  iltered_t  ime_stamps2  )  ; 
f iltered_time_stamps_trans3=transpose  ( f  iltered_time_stamps3 )  ; 

total_t rigger s -neutron s=  [numel  (  f  iltered_t ime -Stamps  1 )  ... 

numel  (  f  iltered_t  ime_stamps2  )  numel  (  f  iltered_time_stamps3  )  ]  ; 

minimum_trigger s_neutrons=min (total_triggers_neutrons ) ; 

display (' Function  Generator  Time  Stamps  Filtered') 


9'2'9-2'2'9'9'9'9'9'9'9-2'2-2'9'9'9'9'9'9'9'2-2'9-9'2'9'9'9'9'9'9'2-2'2'2'9'9'9-9'9'9'2'2-2'2'9'2'9'9'9'9'9'2'9-2'2'2'9-9'9'9'9'9'2'2-2'2- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%  Filtering  neutron  datasets  to  pull  out  the  correct  PSD  value 

2'2'9-2'2-2'2'9'2'2'9-9'2'9'2'2-2'2'9'9'2'9'2-2'9-2'9-9'9'9'9'9'2'9-2'9-2'2-9'2'9'9'2'2'9-2'2-2'9-9'9-2'9'9'2'9-2'2-2'2'9'2-9'2--9'2'9-2'9- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%Only  look  at  neutron  events  for  first  digitizer 
incr2=l ; 

for  i=l : numel (PSD.neutronsl ) 

for  incrl=l: numel (PSD.neutronsl { i} ) 

if  PSD_neutrons 1 { i} ( incr 1 )  >=0 . 1  &&  PSD_neutrons 1 { i} ( incr 1 ) <=0 . 3 
time_stamp_neutronsl  (incr 2)  =timing_neutronsl{i}  (incrl)  ; 
neutronsl ( incr2 ) =PSD_neutrons 1 { i } (incrl) ; 

incr2=incr2+l ; 
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end 


end 

end 

time_stamp_neutronsl_trans=transpose  (time_stamp_neutrons  1 )  ; 
neutronsl_trans=transpose (neutronsl)  ; 

%Only  look  at  neutron  events  for  second  digitizer 
incr2=l ; 

for  i=l : numel (PSD_neutrons2 ) 

for  incrl=l : numel (PSD_neutrons2 { i} ) 

if  PSD_neutrons2 { i} ( incr 1 ) >=0 . 1  &&  PSD_neutrons2 { i} ( incrl )  <=0 . 3 
time_stamp_neutrons2  (incr 2)  =timing_neutrons2{i}  (incrl)  ; 
neutrons2 (incr2) =PSD_neutrons2{i} (incrl) ; 

incr2=incr2+l ; 

end 

end 

end 

time_stamp_neutrons2_trans=transpose  (time_stamp_neutrons2 )  ; 
neutrons2_trans=transpose (neutrons2)  ; 

%Only  look  at  neutron  events  for  third  digitizer 
incr2=l ; 

for  i=l :nurael (PSD_neutrons3) 

for  incrl=l: numel (PSD_neutrons3{i}) 
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if  PSD_neutrons3{ i} ( incrl ) >=0 . 1  &&  PSD_neutrons3{ i} ( incrl ) <=0 . 3 
time_stamp_neutrons3  (incr2)  =timing_neutrons3{  i}  (incrl)  ; 
neutrons 3  ( incr2 ) =PSD_neut rons3 { i } (incrl ) ; 

incr2=incr2+l ; 

end 

end 

end 

time_stamp_neutrons3_trans=transpose  (time_stamp_neutrons3)  ; 
neutrons3_trans=transpose (neutrons3)  ; 


display ('PSD  Values  Obtained') 

9'2'9-9'2'9'9'9'9-9'9'9-9'2-9'2'9'9'9'9'9'9'2-9-2-9'2'9'9-9'9'9'9'2-2'2'2'9'9'9-9'9'9'2'2-2'9'2'2'9'9'9'9'9'9'9-9'2'9'9'9'9-9'9'9'9'2-2'2'9- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%  Compare  frequency  generator  time  stamps  to  neutron  event  time  stamps 
%  Prevent  neutron  data  stamps  to  be  applied  to  incorrect  time  stamps 
%  due  to  clock  turnover 

9'2'2-9'2'9'9-9'9-9'9'2-9'2-2'2'2'9-9'9'9'9'2-2-2-9'2'9'9'2'9'9'9'2-2'2'2'2'9'9-9'9'9'9'2-2'2'9'2'9'9-9'9'9'2'2-9'2'2'9-9'9-9'9'9'9'2-2'2'9- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


t  ime_between_neut  ron_event  s  1  =  0  ; 
old_timeindex_of  .neutron  s  1  =  0; 

for  incr=l : numel (time_stamp_neutronsl ) 

for  i=l :numel (filtered_time_stampsl) 

timing_dif  ference_NaI_neutronsl=time_stamp_neutronsl  (incr)  - 
f iltered_time_stampsl  (i)  ; 

t ime _bet wee n_neutron_e vent s l=abs  ( i-old_timeindex_of -neutrons  1 )  ; 
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if  abs  (t  iming_dif  f  erence_NaI_neutrons  1 )  <= 
concident_acquisition_window 

timeindex_of_neutronsl  (incr)=i; 
old_t  imeindex-of -neutrons l=i ; 
break 

end 


end 

end 

display (' Absolute  Time  Stamps  of  Neutronsl  Determined') 


time_between_neutron_events2=0; 
old_timeindex_of  .neutron  s  2  =  0; 

for  incr=l : numel (time_stamp_neutrons2 ) 

for  i=l : numel ( f iltered_time_stamps2 ) 

timing_dif  f  erence_NaI_neutrons2=time_stamp-neutrons2  (incr)  - 
f iltered_time_stamps2  (i)  ; 

time_between_neutron_event  s2=abs  (i-old_t  imeindex-of -neutrons  2  )  ; 

if  abs  (timing-dif  f  erence_NaI_neutrons2  )  <= 
concident_acquisition_window 

timeindex_of _neutrons2  (incr)=i; 
old_timeindex-Of -neutron  s  2 =i ; 
break 

end 
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end 


end 


display (' Absolute  Time  Stamps  of  Neutrons2  Determined') 


time_between_neutron_events3=0; 
old_timeindex_of  .neutron  s  3  =  0; 

for  incr=l : numel (time_stamp_neutrons3 ) 

for  i=l : numel ( f iltered_time_stamps3 ) 

timing_dif  f erence_NaI_neutrons3=time_stamp_neutrons3  (incr)  - 
f iltered_time_stamps3  (i)  ; 

time_between_neutron_event s3=abs  (i-old_timeindex_of meutrons 3 )  ; 

if  abs  (timing_dif  f  erence_NaI_neutrons3  )  <= 
concident_acquisition_window 

timeindex_of_neutrons3  (incr)  =i; 
old_timeindex_of  _neutrons3=i ; 
break 

end 


end 

end 

display (' Absolute  Time  Stamps  of  Neutrons3  Determined') 

9'2'9-2'9-2'9'9'9'9'9'9'2'9'2'2'2'9'9'9'9'9'9-2'9-2'9-9'9'9'9'9'2'9-2'9-2'2'9'9'9'9'9'2'9-2'9-2'2-9'9'9'9'9'2'9-2'9-2'9-9'9-9'9'9'9'9-2'2'2'9- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 
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Combine  all  neutron  events  time  stamps  into  one  vector 


9'2'9'9'2'9'9'9'9'9'9'9-2'2-9'9'9'9'9'9'9'9'2-9'9-9'2'9'9-9'9'9'9'2-2'2'9'9'9'9-9'9'9'2'2-2'9'9'2'9'9'9'9'9'2'2-2'2'9'9-9'9'9'9'9'2'2-2'2'9'9- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


timeindex_of_neutronsl_trans=transpose  (timeindex_of  .neutron si )  ; 
timeindex.of  .neutron s 2 _t ran s=t ran spose  (timeindex.of  .neutrons 2  )  ; 
timeindex.of  .neutron  s  3  _t  ran  s=t  ran  spose  (timeindex.of  .neutrons 3  )  ; 

total_neutron_indexes=vertcat  (timeindex.of.neutronsl.trans ,  .  .  . 
timeindex.of  .neutron  s  2  _t  ran  s  ,  timeindex.of  .neutron  s  3  _t  ran  s  )  ; 


9'2'2-2'2'2'9-9'2'9'9'2-2'2-2'2'2'2'9'9'9'9'2-2'2-2'2'9'9'9'9'9'9'2-2'2'2-2'9'9-9'9'9'2'2-2'2'2'2'2'2-9'9'9'2'2-2'2-2'9-9'2'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%Save 

%This 

%times 


previous  section  workspace  as  'Material ' .Run.Neutron.Timing 
allows  the  energy  determination  section  to  be  ran  multiple 
without  having  to  rerun  the  nutron  timing  section. 


2'2'2-2'2'2'9'9'9'2'9'2-2'2-2'2'2'2'9'9'9'9'2-2'2-2'2'9'9'9'9'2'2'2'2'2'2'2'2'9'9'2'9'2'2-2'2'2'2'9'9'9'9'9'2'2-2'2'2'9-2'9'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%  clear; 

%  load ( ' Pb3_Run_Neutron_Timing  '  )  ; 
%  clc; 


9'2'9-2'9-2'9-9'9-9'9'9'2'9'2'9-9'9-9'9'9'9'9-2'9-2'2'9'9-9'2'9'9'9-9'2'2'9-9'9-9'9'9'2'9-2'9-2'9-9'9-9'9'9'2'9-2'2-2'9-9'9-9'9'9'2'9-2'9-2' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%  Begin  Nal (Tl)  energy  determination 


9'2'2-9'2'9'9-9'9'9'9'9-2'2-9'2'9'9'9'9'9'9'2-2'9-9'2'9'9'9'9'9'9'2-2'2'2'9'9'9-9'9'9'2'2-2'2'9'2'9'9'9'9'9'2'9-2'2'2'9-9'9'9'9'9'2'2-2'2'9- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


9'2'9'2'9-2'9-9'9-9'9'9-9'9-2'9-2'2-9'9'9'9'9-9'9-9'9-9'9-9'9'9'9'9-2'2'2'9'9'9'9'9'9'2'9-2'9-2'9-9'9-9'2' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%  Load  Bethe  Bloch  Calculation  Data 


2'2'2-2'2'2'9'9'9'2'9'2-2'2-2'2'2'2'9'9'9'9'2-2'2-2'2'9'2'2'9'2'9'2-2'2'2'2'2'9-9'2'9'2'2-2-2'2'2'9'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooo 


load ( 1 BetheBlock . mat 1 ) ; 
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counter  =  1; 
counterl  =  1; 


2'2'2-2'2'2'9'2'9-2'9'2-2'2-2'2'2'2'9'9'9'9'2-2'2-2'2'9'9'2'9'9'9'2-2'2'2'2'2'9-9'9'9'2-2-2'2'2'2'9'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%  Nal  (Tl)  detector  caliberation  data 


9'2'9-2'9-2'2'9'2'2'2-9-2'9-2'2'2'2'9'9'2'9'2-2'9-2'9'9'2'9'9'2'2'9-2'9-2'2'9'2'9'9'9'2'2-2'2-2'9-9'2'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%Used  Cs-137  and  Co-60  peaks  for  calibration 


energy  _NaI_cal  =  [0.662 
MIP  =  [9910,  9459,  8844 


1.173,  1.333,  24.297]; 
11489] ; 


9'2'2-2-2'2'9-9'2'9'9' 

ooooooooooo 


%  15  Jan  % 


9-2'9'2'9-2'2-2'2'9'9' 

ooooooooooo 


channel_NaI_cal { 1 }  = 

[409, 

722, 

814, 

MIP  (1)  ]  ; 

channel-Nal.cal {2 }  = 

[372, 

641, 

728, 

MIP  (2) ]  ; 

channel_NaI_cal{3}  = 

[312, 

542, 

609, 

MIP  (3) ]  ; 

channel_NaI_cal { 4 }  = 

[735, 

1271, 

1439 

i,  MIP  ( 4  )  ] 

pulse_height  =  cell2mat (pulse_heights_NaI ) ; 

pulse_height (pulse_height  >  1.8*10~4|  pulse_height  <  0)  =  NaN; 

[r,c]  =  find (isnan (pulse.height ) ) ; 
pulse_height  (r,  : )  =  0; 

for  i  =  1:4 

Caliberationji}  =  polyf it ( channel_NaI_cal{ i} ,  energy_NaI_cal , 1 ) ; 
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Nal_energy ( : , i )  =  pulse.height ( : , i) *Caliberation{i} ( 1 )  + 

Caliberationji} (2) ; 


end 


2'2'9-2'9-2'2'9'2-9'2-9'2'2'2'2-2'9-9'9'2'9'2'2'2'2'9-2'2'9'9'2'2--9-2'9-2'9-9'2'9'9'9'2'9-2'2-2'9-2'2-9'9'2'2'9-2'9-2'2'9'2-2'9'2'2'2- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%Histogram  of  Landau  Distribution  (Energy  Deposited  in  each  Nal(Tl) 


9'2'2-2'2'2'9-2'9-9'9'2-2'2-2'2'2'2-9'9'9'9'2-2'2-2'2'9'9-9'9'9'9'2-2'2'2'2'9'9-9'9'9'2'2-2'2'2'2'9'9'9'9'9'2'2-2'2'2-9'2'9-9'9'9'2'2- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


figure 

for  i  =  1:4 


subplot (2, 2, i) 

hist (Nal_energy ( :  ,  i )  ,  [0:1:60]) 

%  ylim( [0, 1000] ) 
xlim ( [0, 50] ) 

title ([ ’Nal  (Tl)  '  num2str(i)],  'FontSize',  30.) 
xlabel ( ' Energy  Deposited ' ,  ' FontSize ' , 30  .  ) 
ylabel ( ' Counts ' ,  1 FontSize  '  , 30  . ) 
set (gca,  'FontSize',  20) 

end 


9'2-2'9'2'2'9'9'9'9'9'2'2-2'2'9'9'9'9'9'9'9'2'2-2'2'9'9'9'2'9'9'9'9'2'9'2'9'9'9'9-9'9'9'2'2-9'9'9'9'9'9'9'9'2-9'2-9'9'2'9-9'9-9'9'2-2-2-9'2'9'9'9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%  Determine  Limits 


to  Muon  Energy 


From  Deposited  Energy 


9'2-2'2'9'9'9'9'9'9'9'2'2-2'2'2'2'9'9-9'9'9'9'2-2'2'9'9-9'9-9'9'9'2'2-9'2'2'9-2'9-9'9'9'9'2-9'2'9'9'9'9'9'9'2-2'2-9'2'2'9'9'9'9'9'2-9-2-2'2'9'9-9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


lowest_energy=24 . 297  *  0.96; 


%  MIP  24.297 


upper_energy_detectorsl23=31 . 7466  *  1.12; 


%highest  value 


for  5  GeV 
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muon=31 .7466 


upper_energy_detector4=31 . 7466  *  1.12; 


for  m=l : minimum.triggers.neutrons 


%  Set  the  energy  deposited  by  the  muon 


energy l_interp= 
energy2_interp= 
energy3_interp= 
energy4_interp= 


Nal.energy  (m,  1) ; 
Na I .energy  (m,  2 ) ; 
Nal.energy  (m,  3) ; 
Nal.energy  (m,  4 ) ; 


if  energy l_interp<=upper_energy_detector s  12 3  &&  energy2_interp<= 

upper_energy_detectorsl23  .  .  . 

&&  energy3_interp<=upper_energy_detectorsl23  && 
energy4_interp<=upper_energy_detector4  .  .  . 

&&  energyl_interp>=lowest_energy  &&  energy2_interp>= 
lowest.energy . . . 

&&  energy3_interp>=lowest_energy  &&  energy4_interp>= 
lowest.energy 


9'2'9'2'9-2'9-2'9'9'9-9'9'9'2'9'2'9-2'9-9'9-9'9'9'2'9-2'9-2'9-9'9-9'9'9'9'9'2'9-2'9-9'9-9'9'9'2'9'2'9-2'9-9'2'9'9'9-2'9-2'2-2'9-9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%Count  all  neutrons  for  that  muon  event. 


o_ 

o 


9'9'9'2'2-2'2'2'2'9'9'9'9'9'2'2-2'2'2'9'9'9'9'9'9'2'2-2'2'2'9'9'9-9'9'9'2'2-2'2'2'9-9'9-9'9'2-2-2-2'2'2'9'9'9'9'9'2-9-2-2'2'9'9-9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 
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neutrons=0 ; 


for  j=l : numel (total_neutron_indexes ) 
if  m==total_neutron_indexes ( j ) 
neutrons=neutrons+l ; 


end 

end 

Find  Chi~2 

chi 2= ( (dElBB_NaI-energyl_interp) . "2 . / dEIBB.Nal ) + . . . 

(  (dE2BB_NaI-energy2_interp)  . ~2./dE2BB_NaI)+.  .  . 

( (dE3BB_NaI-energy3_interp) . "2 . / dE3BB_NaI ) + . . . 

( (dE4BB_NaI-energy4_interp) .  ~2  .  / dE4BB_NaI ) ; 

chi2_list=transpose (vertcat (chi2,  KEnergy_exit_BB ) ) ; 


Determine  the  Lowest  Chi'2 
f  chi2~=0 

[chi2_min, index] =min (chi2_list ( : , 1) ) ; 
chi_list (counterl ) =chi2_min; 

KEnergy_exit=chi2_list (index, 2) ; 


%Determine  the  exiting  energy 
energy4_exiting_list (counterl) =KEnergy_exit ; 
energy4_exiting_hist= [energy4_exiting_hist, KEnergy_exit 
]  ; 
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neutrons_list (counterl) =neutrons; 


counter l=counter 1+1 ; 


end 


end 

%View  Progress  Through  Files  and  clear  chi2  for  next  file 

count er=counter+l ; 
disp (counter) 
chi2=0 ; 
incr 1=1 ; 
chi2_list= [ ] ; 


end 


energy4_exiting_list_trans=transpose  (energy4_exiting_list)  ; 
neutrons_list_trans=transpose (neutrons_list )  ; 

NeutronDataFile= ' NeutronDataFile . txt '  ; 

M=  [energy4_exiting_list_trans  neutrons_list_trans  ]; 

dlmwrite (NeutronDataFile, M,  ' delimiter '  ,  ' \t '  ,  'precision',  ' %15 . 1  Of ' ) 


9'2'9-2'9-2'9-9'9-9'9'9-9'9-2'9'2'9-9'2'9'9'9-2'9'2'9-9'9-9'9'9'2'9-9'9-2'9'9'9-9'9'9'2'9- 

ooooooooooooooooooooooooooooooooooooooooooooo 


9'2'9'9'2'2'9'9'9-9'9'9-2'2-9'2'9'9'9'9'9'9'2-2'9'9'9'9'9'9'9'9'9'2-2'2'2-9'9'9-9'9'9'2'2- 

ooooooooooooooooooooooooooooooooooooooooooooo 


%ANALYSIS  portion  of  code 
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Written  by  2Lt  Logan  Brandt 


clear; 

clc; 

close  all; 


%%%Index  m  =  1  2 

Material  =  cellstr ({ ' BackgroundF ' ,  'FeF 
Color  =  cellstr ({' k !  ,  'b',  ' r',  ' g ' } )  ; 

Mater ial_Name  =  cellstr ({' Background ' , 


3  4 

,  ' PblF 1 ,  ' Pb3F  '  })  ; 

15  cm  Iron',  '5  cm  Lead', 


Lead ' } ) ; 


'15  cm 


bins  =  0:40; 

centers  =  { [ 0 : 1 00 : 5000 ]  ,  [ 0 : 1  :  30 ] } ; 


for  m  =  1 : length (Material ) 


9'9'9-2'9-2'9-9'9-9'9'9'2'9'2'9-2'9-9'9'9'2'9-2'9'2'2'9'9-9'2'9'9'9-9'9'2'9-9'9-9'9'9'2'9-2'9-9'9-9'2-9'9'9'2-9-2'9-2'9-9'9-9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%Load  Data 


9'2'2'2'2'2'9'9'9'2'9'2-2'2-2'2'2'2'9'9'9'9'2-2'2-2'2'9'2'2'9'9'9'2-2'2'2'2'2'9-9'9'9'2'2-2'2'2'2'9'9'9'9'9'2'2-2'2-2'9'2'9'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

Data{m}  =  load ( [Materialjm}  '_NeutronDataFile.txt']); 
Muon_Energy{m}  =  Datajm} ( : , 1 ) ; 

Neutron_Count{m}  =  Datajm} ( :  ,  2 ) ; 

Total_Neutrons{m}  =  sum (Neutron_Count{m} ) ; 

Total  JYIuonsjm}  =  length  (Muon_Energy{m}  )  ; 
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Muon  Energy  Distribution 


2-2'2'2-2'2-9'2-9'9--9-2'9-2'9-2'2-2'9'2'9--9-2'9-2'2-9'2-2'9'9'9'9-2'9-2'9-9'2-9'9--9'2'9-2'9-2'2'9'2-9'9'9'2'9'2'9-2'2-9'2'9'9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


Energy_Hist{m}  =  hist (Muon_Energy{m} ,  centersjl}); 
Mean_Muon_Energy{m}  =  mean  (Muon_Energy {m} )  ; 
Mode_Muon_Energy{ra}  =  mode  (Muon.Energy {m} )  ; 


2'2-9'2'9'9'9'9-9'9'2-9'2-2'9'9'9-9'9'9'9'2-2'2-9'2'9'9-9'9'9'9'2-9'2'2'2'9'9-9'9'9'9'2-2'9'9'2'9'9-9'9'9'2'9-2-2'2'9'9'9-9'9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

Neutron  Muon  Energy  Dependence 

2'2-2'2'2'2-9'2-9'9'2-2-2-2'2'2'9-2'9'9'9'2-2'2-2'2'9'9-9'9'9'9'2-2'2'2'2'9'9'9'9'9'2'2-2'2'2'2'9'9-9'9'9-2'2-2'2'2'9-9'9-9'9'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

Hist2d{m}  =  hist3 (Data{m} ,  'Ctrs',  centers,  'FaceColor',  Colorjm}) 

for  r  =  1 : length (centersjl}) 

Average_Neutrons_per_Energy{m} ( r )  =  0; 

for  c  =  1 : length (centers{2} ) 

Average_Neutrons_per_Energy{m}  (r)  = 

Average_Neutrons_per_Energy{m}  ( r )  +  Hist2d{m} ( r, c) * (c-1 ) 


end 

Neutrons_per_Energy{m}  (  r )  =  Average  _Neutrons_per_Energy{m}  (  r )  ; 
Muons_per_Energy{m} ( r )  =  sum (Hist2d{m} ( r , : ) ) ; 

Average  _Neutrons_per_Energy{m}  (  r)  =  Neutrons_per_Energy{m}  (  r)  / 
Muons_per_Energy{m} (r) ; 

Average_Neutrons_per_Energy_Error  {m}  (  r )  =  sqrt((sqrt( 

Neutrons_per_Energy {m}  (r)  )  /Neutrons_per_Energy{m}  (r)  )  ~2  .  .  . 
+  (sqrt (Muons_per_Energy{m} (r) ) /Muons_per_Energy {m} (r) ) ~2) * 
Average_Neutrons_per_Energy{m}  (r)  ; 


end 


for  c  =  1 : length (centers{2}) 

Muons_per_Neutron_Count{m} (c)  =  sum (Hist2d{m} ( : , c) )  ; 
Energy_Spectrum_per_Neutron_Count{m}  (c,  :  )  =  Hist2d{m} ( : , c) ; 


end 


9'2'2-2'2'2'9'9'9-9'9'2-2'2-2'2'2'2'9'9'9'9'2-2'2-2'2'9'2'9'9'9'9'2-2'2'2'2'9'2'9'9'9'2'2-2'2'2'2'9'9'9'9'9'2-2-2'2'2'9-2'9-9'9'2' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%Bin  Neutrons 


9'2'2-2'2'2'9-9'2-9'9'2-2'2-2'2'2'2'9'9'9'9'2-2'2-2'2'2'9-9'9'9'9'2-9'2'2'2'2'9-9'9'9'2'2-2'2'2'2'9'9-9'9'9'2'2-2'2-2'2-9'2'9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


for  i  =  0:40 

Histjm} ( i+1 )  =  sum (Neutron.Count {m}  ==i); 

end 


Histjm} 


transpose (Histjm} )  ; 


9'2'2-9'2'9'9'9'9-9'9'9-2'2-2'2'9'9'9'9'9'9'2-2'2-2'2'9'9-9'9'9'9'2-9'2'9'9'9'9'9'9'9'9'2-2'2'9'2'9'9'9'9'9'2'2-9-2'2'9'9'9-9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%Remove  Background 

9'2'9-2'9-2'2-9'9-9'9'9-9'9-2'9-2'9-9'9'9'2'9-2-9'2'2'9'9-9'2'9'9'9-9'9'2'9-9'9'9'9'9'9'9-2'9-2'2'9'9-9'2'9'2'9-2'9'2'9-9'9-9'9'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

Hist_above_Background{m}  =  Hist{m}  -  Histjl}; 

Hist_above_Background{m}  (Hist_above_Background{m}  (  :  )  <  0)  =0; 
Muons_above_Background{m}  =  sum (Hist_above_Background{m}  )  ; 
Neutrons_above_Background{m}  =  0; 
for  i  =  0:40 

Neutrons_above_Background{m}  =  Neutrons_above_Background{m}  + 
Hist_above_Background{m} (i+1) *i; 

end 
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%Find  Results  and  error  to  produce  specific  number  of  neutrons 


2'2'9-2'9-2'2-2'2'9'9-9-2'9-2'2'2'9-9'2'9'2'2'2'9-2'9'9'2-9'9'2'2'9-2'9-2'2-9'2'9'2'2'2'9-2'2-2'9-9'2-2'9'2'2'9-2'9-2'2-9'2'9'2'2'2' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%Probability 

Probabilityjm}  =  Hist_above_Background{m}  .  /Muons_above_Background{m 

}; 

Prob_Error{m}=sqrt ( (sqrt (Hist_above_Background{m} )  / 
Hist_above_Background{m} ) ~2  ... 

+  (sqrt  (Muons_above_Background{m}  )  /Muons_above_Background{m} )  "2) 
*Probability{m} ; 

%Average  Neutrons  produced  found  from  total  numbers 
Average_Neutrons{m}  =  Total_Neutrons{m}/Total_Muons{m}; 
Average_Neutrons_Error {m}  =  sqrt ( (sqrt (Total_Neutrons{m}  ) / 
Total.Neutrons {m} ) "2  ... 

+  (sqrt (Total_Muons {m} ) /Total-Muons {m} ) ~2)  *Average_Neutrons{m} ; 

%Expectation  value  for  the  neutron  count  per  muon  found  with 
background 
%removed 

Neutrons_per_Muon{m}  =  Neutrons_above_Background{m}  / 
Muons_above_Background{m} ; 

Neutrons_per jyiuon_Error{m}  =  sqrt  (  (sqrt  (Neutrons_above_Background{m 
})  /Neutrons_above_Background{m} )  "2  ... 

+  (sqrt  (Muons_above_Background{m}  )  /Muons_above_Background{m} )  "2) 
*Neutrons_per_Muon{m} ; 


end 
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o,  o. 
o  o 


for  m  =  1 : length (Material ) 

for  i  =  1:  length  (Average_Neutrons_per_Energy{m} ) 
hO(m, i)  =  ttest (Average_Neutrons_per_Energy{m} (i) - 
Average_Neutrons_per_Energy{l} (i) ) ; 

end 

MeV_Avg{m}  =  ( Average_Neutrons_per_Energy{m}  ( 2  )  + 
Average_Neutrons_per_Energy{m}  (3)  )  /2; 

GeV_Avg{m}  =  (Average_Neutrons_per_Energy{m}  (20)  + 

Average_Neutrons_per_Energy{m}  (21)  )  / 2; 

MeV_Avg_Er  ror  {m}  =  (Average_Neutrons_per_Energy_Error  {m}  ( 2  )  + 
Average  _Neutrons_per_Energy_Error{m}  (3)  )  /2; 

GeV_Avg_Error{m}  =  (Average_Neutrons_per_Energy_Error{m}  (20)  + 
Average_Neutrons_per_Energy_Error{m}  (21)  )  /2; 

end 

for  m  =  2:  length (Material ) 

MeV_Avg{m}  =  MeV_Avg{m}  -  MeV_Avg{l}; 

GeV_Avg{m}  =  GeV_Avg{m}  -  GeV_Avg{l}; 


end 


9'2'2-2'2'2'9-2'9-9'9'2-2'2-2'2'9'9'9'9'9'9'2-2'2-9'2'9'9-9-9'9'9'2-2'2'2'9'9'9- 

oooooooooooooooooooooooooooooooooooooooo 


%T-Tests  comparing  data  sets 


9'2'2-2'2'2'9-2'9-9'9'2-2'2-2'2'2'2'2'9'9'9'2-2-2-2'2'9'2'2'9'9'9'2-2'2'2'2'9'2- 

oooooooooooooooooooooooooooooooooooooooo 
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[hl,pPblPb3]  =  ttest2 (Hist_above_Background{3} ,  Hist_above_Background 
{4}) 

[h2,pFePb3]  =  ttest2 (Hist_above_Background{2} ,  Hist_above_Background{4} ) 
[h3,pFePbl]  =  ttest2 (Hist_above_Background{2} ,  Hist_above_Background{3}) 

[hl,pBkgFe]  =  ttest2 (Hist_above_Background{ 1} ,  Hist_above_Background{2}) 
[h2,pBkgPbl]  =  ttest2 (Hist.above.Backgroundjl} ,  Hist.above.Background 

{3}) 

[h3,pBkgPb3]  =  ttest2 (Hist.above.Backgroundj 1} ,  Hist.above.Background 
{4}) 

o,  o. 
o  o 

[hl,pPblPb3]  =  ttest2  (Average_Neutrons_per_Energy{3}  , 

Aver age_Neutr on s_per .Energy  {  4  }  ) 

[h2,pFePb3]  =  ttest2  (Average_Neutrons_per_Energy{2}  , 

Ave r age _Neutr on s_per .Energy  {  4  }  ) 

[h3,pFePbl]  =  ttest2 (Average_Neutrons_per_Energy{2}  , 

Average  .Neutrons  .per  .Energy  {3}  ) 

[hl,pBkgFe]  =  ttest2 (Average.Neutrons.per.Energyjl}  , 

Average.Neutrons .per .Energy  {2  }  ) 

[h2,pBkgPbl]  =  ttest2  (Average_Neutrons_per_Energy{l}  , 
Average_Neutrons_per_Energy  {3} ) 

[h3,pBkgPb3]  =  ttest2  (Average.Neutrons.per.Energyjl}  , 
Average_Neutrons_per_Energy  {  4  }  ) 


9'2'2-2'2'2'9-2'2'9'9'2-2'2-2'2'2'2'9'2'9'9'2-2'2-2'2'9'9'2'9'9'9'2-2'2'2'2'9'2- 

oooooooooooooooooooooooooooooooooooooooo 


%Plot  raw  neutron  counts  per  muon  event 


9'2'2-2'2'2'9-2'9-9'9'2-2'2-2'2'2'9-9'9'9'9'2-2'2-2'2'9'9'9'9'9'9'2-2'2'2-2'9'2- 

oooooooooooooooooooooooooooooooooooooooo 


figure ( 1 ) 
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groupplot=  horzcat (Histjl : length (Material) }) ; 

Neutron_Histogram  =  bar (bins,  [groupplot],  'grouped'); 
xlim ( [-1,21] ) 

for  m  =  1 : length (Material ) 

set (Neutron_Histogram (m) ,  'FaceColor',  Color{m},  'EdgeColor',  Color 
{m})  ; 

end 

legend (Mater ial_Name ) 

xlabel (' Neutron  Counts  per  Muon '  ,  ' FontSize '  ,  30 . ) 
y label ( ' Muon  Events ', 'Font Size', 30.) 
set (gca,  'FontSize',  20) 


9'2'2-2'2'2'9-2'9'9'2'2-2'2-2'2'2'2'9'9'9'9'2-2'2-2'2'9'2'9'9'9'9'2-2'2'2'2'9'2'9'9'9'2- 

oooooooooooooooooooooooooooooooooooooooooooo 

%Plot  Neutron  Counts  with  Background  Removed 

2'2'9-2'9'2'9-9'9'9'9'9-9'9-2'9-2'9-9'9'9'9'9-9-9-2'2'9'9-9'9'9'9'9-9'2'2'9-9'9-9'9'9'2' 

oooooooooooooooooooooooooooooooooooooooooooo 

figure ( 2 ) 

groupplot=  horzcat (  Hist_above_Background{2 : length (Material) }) ; 
Neutron_Histogram_wo_Background  =  bar (bins,  [groupplot],  'grouped'); 
xlim ( [-1,21] ) 

for  m  =  2 : length (Material ) 

set (Neutron_Histogram_wo_Background (m-1 ) ,  'FaceColor',  Colorjm},  ' 
EdgeColor',  Colorjm}) ; 

end 

legend (Material_Name{2 : length (Material-Name ) }  ) 
xlabel (' Neutron  Counts  per  Muon ',' FontSize ', 30  .  ) 
ylabel('Muon  Events',  ' FontSize ', 30  .  ) 
set (gca,  'FontSize',  20) 
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%Plot  Probability  Distribution  with  Error  Bars 

2'2'9-2'9-2'2'9'2'9'9-9-2'9-2'2'2'2-9'9'2'9'2-2'9-2'9-9'9'9'9'2'9'9-2'9-2'2'9'2'9'2--2'2'9-2' 

oooooooooooooooooooooooooooooooooooooooooooooo 

figure (3) 

for  m  =  2 : length (Material ) 

%  Poisson_plot{m}  =  plot (bins , poisspdf (bins ,  Neutrons_per_Muon{m 

}) ,  Color{m}) ; 

Prob_Plot{m}  =  errorbar (bins ,  Probability{m} ,  Prob_Error {m}  , 
Colorjm},  ' MarkerSize ' ,  10); 
hold  all 


end 

poissfit (Probability{m} ) 
xlim ( [-1,21] ) 
ylim( [0, .25] ) 

legend (Material_Name{2 : length (Material-Name ) }  ) 
xlabel (' Neutrons  Produced  per  Muon ' , ' FontSize ' , 30 . ) 
y label ( ' Probability  1 ,  ' FontSize ' , 30 . ) 
set (gca,  'FontSize',  20) 


9'2'2-2'2'2'9-9'2-9'9'2-2'2-2'2'2'2'9'9'9'9'2-2'2-2'2'9'2'9'9'9'9'2-2'2'2'2'9'2- 

oooooooooooooooooooooooooooooooooooooooo 


%Plot  Neutron  counts  against  Muon  Energy 


9'2'9-2'9-2'9-9'9-9'9'9'2'9'2'9-2'9-9'9'9'9'9'2'9'9'2'9'9-9'2'9'9'9-2'2'2'9-9'9- 

oooooooooooooooooooooooooooooooooooooooo 


figure (5) 


for  m  =  1 : length (Material ) 
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subplot (2 , 2 , m) 

scatter (Muon_Energy{m},  Neutron_Count{m},  ’  .  ' MarkerEdgeColor Color 


{m} . . . 

, ' MarkerFaceColor '  ,  Color{m},'SizeData',250); 
title (Mater ial_Name{m} ,  'FontSize',  30.) 
xlim ( [0, 3000] ) 

xlabel ( 1  Muon  Energy  [MeV]  ' ,  ' FontSize ' , 30  .  ) 
ylim( [0,30] ) 

ylabel ( ' Neutron  Counts ' ,  ' FontSize ' , 30  .  ) 
hold  on 

set (gca,  'FontSize',  20) 

end 


9'2'2-2'2'2'9-2'9-9'9'2-2'2-2'2'2'2'9'9'9'9'2-2'2-2'2'9'2'9'9'9'9'2-2'2'2'2'9'9-9'9'9'2'2-2'2'2-2'9'9-9'9'9'2'2- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%Neutron  Yield  for  each  100  MeV  incident  muon  energy  bin 


9'2'9'2'9-2'9-9'9-9'9'9'2'9'2'9-2'9-9'9'9'9'9-2'9-2'2'9'2-9'9'9'9'9-2'2'2'9-9'9-9'9'9'2'9-2'9-2'9-9'9-9'2'9'2'9- 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


figure ( 6) 

for  m  =  1 : length (Material ) 

subplot (2 , 2 , m) 

bar ( centers{ 1} ,  Average_Neutrons_per_Energy{m}  ,  'FaceColor' 

m} ,  .  .  . 

'EdgeColor',  Colorjm}); 
ylim( [0, 10] ) 
xlim ( [0, 3000] ) 

title (Material_Name{m} ,  'FontSize',  30.) 
xlabel ('Muon  Energy  [MeV] ',' FontSize ', 30 . ) 
ylabel (' Average  Neutron  Counts ',' FontSize ', 20  .  ) 


Colorj 
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set (gca,  'FontSize',  20) 


end 


2'2'2-2'2'2'9'9'9-2'9'2-2'2-2'2'9-2'9'9'9'9'2-2'2-2'2'9-9'2'9'2'9'2- 

oooooooooooooooooooooooooooooooooo 


%Incident  Muon  Energy  Distribution 


9'2'9-2'9'2'2'9'2'2'9-9-2'2-2'2'2'2-9'2'2'9'2'2'9-2'9-9'2'9'9'2'2--9- 

oooooooooooooooooooooooooooooooooo 


figure ( 7 ) 

for  m  =  1 : length (Material ) 

subplot (2, 2,  m) 

bar ( centers{ 1} ,  Energy_Hist{m} ,  'FaceColor',  Color{m}, 
Color{m} ) ; 
ylim ( [0,4000] ) 
xlim ( [0, 3000] ) 

title (Mater ial_Name{m} ,  'FontSize',  30.) 
xlabel('Muon  Energy  [MeV] ',' FontSize ', 30 . ) 
y label ( ' Muon  Counts ' ,  ' FontSize ' , 30  . ) 
set (gca,  'FontSize',  20) 


end 


2'2'2-2'2'2'9'9'2-9'9'2-2'2-2'2'2'2'9'9'9'9'2-2'2-2'2'9'9'2'9'2'2'2-2'2'2'2-2'9-9'9'9'2'2-2' 

oooooooooooooooooooooooooooooooooooooooooooooo 


%Plot  Average  Neutron  Yield  with  Error 


2'2'9-2'9-2'2-9'2'2'9-9-2'9-2'2'2'9'9'2'2'9'9'2'9-2'9-9'2-9'9'2'2'9-2'9-2'2-9'2'9'2--2'2'2-2' 

oooooooooooooooooooooooooooooooooooooooooooooo 


1 EdgeColor ' , 
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figure (8) 

for  m  =  1 : length (Material ) 

Avg_Energy_Plot{m}  =  errorbar (centersjl} , 
Average_Neutrons_per_Energy{m} , . . . 

Average_Neutrons_per_Energy_Error{m} ,  Colorjm},  ' MarkerSize ' 
,  10); 
hold  all 

end 

xlim ( [0, 3000] ) 
ylim( [0, 10] ) 

legend (Mater ial_Name{ 1 : length (Material-Name ) }  ) 
xlabel (' Incident  Muon  Energy  [MeV] ' , ' FontSize ' , 30 . ) 
ylabel (' Average  Neutron  Counts FontSize 30  .  ) 
set (gca,  'FontSize',  20) 


9'2'9'2'9-2'9-9'9-9'9'9'2'9'2'9-2'9-9'9'9'2'9-9'9'2'9-9'9-9'9'9'9'9-2'2'2'9-9'9-9'9'9'2'9-2'9-2'2-9'9-9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooo 

%Plot  Muon  Energy  Spectrum  for  a  Given  Neutron  Yield 

9'2'2-2'2'2'9-2'9-2'9'2-2'2-2'2'2'2'9'9'9'9'2-2'2-2'2'9'2-2'9'9'9'2-2'2'2'2'2'9-9'9'9'2-2-2'2'2'2'9'9'2' 

oooooooooooooooooooooooooooooooooooooooooooooooooooo 

Neutron_Count_Analyzed  =  4; 

%  figure ( 9 ) 

for  neutrons  =  1:  Neutron_Count_Analyzed 
subplot (2,2, neutrons ) 
for  m  =  1 : length (Material ) 

Energy_Spectrum_per-Neutron-Count_Plot{m}  =  plot (centersjl} , 

Energy-Spectrum_per-Neutron_Count{m}  (neutrons,  :  )  ,  Color{ 
m} ,  'MarkerSize',  10); 
hold  all 

end 

xlim ( [0, 3000] ) 


98 


legend (Mater ial_Name{ 1 : length (Material_Name ) } ) 
title ([' Neutron  Yield  of  '  num2str (neutrons-1) 
FontSize '  ,  30 . ) 

xlabel (' Incident  Muon  Energy  [MeV] ', 'FontSize 
ylabel('Muon  Counts’,  ' FontSize 30  . ) 
set (gca,  'FontSize',  20) 
end 


'  Neutrons '  ]  , 

,25.  ) 
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